Merge pull request #35 from lcmd-epfl/sad_fix

Fix SPAHM[SAD] for open-shell systems #32

qstack/spahm/guesses.py

@@ -55,7 +55,13 @@ def SAD(mol, func):
   mf = pyscf.dft.RKS(mol)
   mf.xc = func
   vhf = mf.get_veff(dm=dm)
-  fock = hc + vhf
+  if vhf.ndim == 2:
+      fock = hc + vhf
+  else:
+      fock = hc + vhf[0]
+      if not numpy.array_equal(vhf[0], vhf[1]):
+        msg = f'The effective potential ({func}) return different alpha and beta matrix components from atomicHF DM'
+        warnings.warn(msg)
   return fock
 
 def SAP(mol, *_):

tests/test_global.py

@@ -8,7 +8,8 @@ def test_avg_kernel():
 
     path = os.path.dirname(os.path.realpath(__file__))
     X_dir = os.path.join(path, 'data/SPAHM_a_H2O/')
-    mols = [np.load(os.path.join(X_dir, f), allow_pickle=True) for f in os.listdir(X_dir) if os.path.isfile(os.path.join(X_dir,f))]
+    mollist = [os.path.join(X_dir, s) for s in ['X_H2O.npy', 'X_rotated_H2O.npy', 'X_H2O_dist.npy']]
+    mols = [np.load(f, allow_pickle=True) for f in mollist]
     K = kernel.kernel(mols, akernel='L', gkernel='avg', sigma=1.0)
 
     true_K = np.array(  [[1.       ,  1.       ,  0.79179528], \
@@ -22,7 +23,8 @@ def test_avg_kernel():
 def test_rem_kernel():
     path = os.path.dirname(os.path.realpath(__file__))
     X_dir = os.path.join(path, 'data/SPAHM_a_H2O/')
-    mols = [np.load(os.path.join(X_dir, f), allow_pickle=True) for f in os.listdir(X_dir) if os.path.isfile(os.path.join(X_dir,f))]
+    mollist = [os.path.join(X_dir, s) for s in ['X_H2O.npy', 'X_rotated_H2O.npy', 'X_H2O_dist.npy']]
+    mols = [np.load(f, allow_pickle=True) for f in mollist]
     K = kernel.kernel(mols, akernel='L', gkernel='rem', sigma=1.0, gdict={'alpha':1.0, 'normalize':1})
 
     true_K = np.array(  [[1.      , 0.6528238, 1.      ], \

tests/test_spahm_a.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+
+import os
+import numpy as np
+from qstack import compound
+from qstack.spahm.rho import atom
+
+
+def test_water():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'minao', charge=0, spin=None)
+
+    X = atom.get_repr(mol, ["H", "O"], 0, None, dm=None,
+                      guess='LB', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+
+    X_true = np.load(path+'/data/SPAHM_a_H2O/X_H2O.npy', allow_pickle=True)
+    assert(X.shape == X_true.shape)
+    for a, a_true in zip(X, X_true):
+        assert(a[0] == a_true[0])                        # atom type
+        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
+
+def test_water_SAD_guess_open_shell():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=1, spin=1) ## test breaks when effective open-shell caluclation is needed
+
+    Xsad = atom.get_repr(mol, ["H", "O"], 1, 1, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+    Xtrue = np.load(path+'/data/SPAHM_a_H2O/X_H2O-RC_SAD.npy', allow_pickle=True)
+    assert(Xsad.shape == Xtrue.shape)
+    for a, a_true in zip(Xsad, Xtrue):
+        assert(a[0] == a_true[0])                        # atom type
+        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
+
+def test_water_SAD_guess_close_shell():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0) ## test breaks when effective open-shell caluclation is needed
+
+    Xsad = atom.get_repr(mol, ["H", "O"], 0, None, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+    Xtrue = np.load(path+'/data/SPAHM_a_H2O/X_H2O_SAD.npy', allow_pickle=True)
+    assert(Xsad.shape == Xtrue.shape)
+    for a, a_true in zip(Xsad, Xtrue):
+        assert(a[0] == a_true[0])                        # atom type
+        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
+
+
+
+if __name__ == '__main__':
+    test_water()
+    test_water_SAD_guess_close_shell()
+    test_water_SAD_guess_open_shell()