small-fixes

qstack/spahm/rho/bond.py

@@ -10,7 +10,7 @@
 
 def bond(mols, dms,
          bpath=defaults.bpath, cutoff=defaults.cutoff, omods=defaults.omod,
-         spin=None, elements=None, only_m0=False, zeros=False, printlevel=0,
+         elements=None, only_m0=False, zeros=False, printlevel=0,
          pairfile=None, dump_and_exit=False, same_basis=False, only_z=[]):
     """ Computes SPAHM-b representations for a set of molecules.
 
@@ -20,7 +20,6 @@ def bond(mols, dms,
         - bpath (str): path to the directory containing bond-optimized basis-functions (.bas)
         - cutoff (float): the cutoff distance (angstrom) between atoms to be considered as bond
         - omods (list of str): the selected mode for open-shell computations
-        - spin (list of int): list of spins for each molecule
         - elements (list of str): list of all elements present in the set of molecules
         - only_m0 (bool): use only basis functions with `m=0`
         - zeros (bool): add zeros features for non-existing bond pairs
@@ -40,8 +39,6 @@ def bond(mols, dms,
     elements, mybasis, qqs0, qqs4q, idx, M = dmbb.read_basis_wrapper(mols, bpath, only_m0, printlevel,
                                                                      elements=elements, cutoff=cutoff,
                                                                      pairfile=pairfile, dump_and_exit=dump_and_exit, same_basis=same_basis)
-    if np.array(spin==None, ndmin=1).all():
-        omods = [None]
     qqs = qqs0 if zeros else qqs4q
     maxlen = max([dmbb.bonds_dict_init(qqs[q0], M)[1] for q0 in elements])
     if len(only_z) > 0:
@@ -116,10 +113,13 @@ def get_repr(mols, xyzlist, guess,  xc=defaults.xc, spin=None, readdm=None,
         all_atoms   = np.array([z for mol in mols for z in mol.elements if z in only_z], ndmin=2)
     else:
         all_atoms   = np.array([mol.elements for mol in mols])
+    spin = np.array(spin) ## a bit dirty but couldn't find a better way to ensure Iterable type!
+    if (spin == None).all():
+        omods = [None]
 
     allvec  = bond(mols, dms, bpath, cutoff, omods,
-                   spin=spin, elements=elements,
-                   only_m0=only_m0, zeros=zeros, printlevel=printlevel,
+                   elements=elements, only_m0=only_m0,
+                   zeros=zeros, printlevel=printlevel,
                    pairfile=pairfile, dump_and_exit=dump_and_exit, same_basis=same_basis, only_z=only_z)
     maxlen=allvec.shape[-1]
     natm = allvec.shape[-2]
@@ -180,8 +180,8 @@ def main():
 
     if args.filename.endswith('xyz'):
         xyzlist = [args.filename]
-        charge  = [int(args.charge) if args.charge is not None else 0]
-        spin    = [int(args.spin)   if args.spin   is not None else None]
+        charge  = np.array([int(args.charge) if args.charge is not None else 0])
+        spin    = np.array([int(args.spin)   if args.spin   is not None else None])
     else:
         xyzlistfile = args.filename
         xyzlist = utils.get_xyzlist(xyzlistfile)
@@ -195,7 +195,10 @@ def main():
                       elements=args.elements, only_m0=args.only_m0, zeros=args.zeros, split=args.split, only_z=args.only_z)
     if args.print > 0: print(reps.shape)
     if args.merge:
-        np.save(args.name_out+'_'+'_'.join(args.omod), reps)
+        if (spin == None).all():
+            np.save(args.name_out, reps)
+        else:
+            np.save(args.name_out+'_'+'_'.join(args.omod), reps)
     else:
         for vec, omod in zip(reps, args.omod):
             np.save(args.name_out+'_'+omod, vec)

qstack/spahm/rho/utils.py

@@ -53,14 +53,14 @@ def load_mols(xyzlist, charge, spin, basis, printlevel=0, units='ANG', ecp=None,
     return mols
 
 
-def mols_guess(mols, xyzlist, guess, xc=defaults.xc, spin=None, readdm=False, printlevel=0):
+def mols_guess(mols, xyzlist, guess, xc=defaults.xc, spin=[None], readdm=False, printlevel=0):
     dms = []
     guess = guesses.get_guess(guess)
-    for xyzfile, mol in zip(xyzlist, mols):
+    for xyzfile, mol, sp in zip(xyzlist, mols, spin):
         if printlevel>0: print(xyzfile, flush=True)
         if not readdm:
             e, v = spahm.get_guess_orbitals(mol, guess, xc=xc)
-            dm   = guesses.get_dm(v, mol.nelec, mol.spin if spin is not None else None)
+            dm   = guesses.get_dm(v, mol.nelec, mol.spin if sp is not None else None)
         else:
             dm = np.load(readdm+'/'+os.path.basename(xyzfile)+'.npy')
             if spin and dm.ndim==2:

tests/test_spahm_b.py

@@ -9,19 +9,31 @@ def test_water():
     mols = utils.load_mols([xyz_in], [0], [0], 'minao')
     dms = utils.mols_guess(mols, [xyz_in], 'LB', spin=[0])
     X = bond.get_repr(mols, [xyz_in], 'LB', spin=[0], with_symbols=False, same_basis=False)
-    #X = np.hstack(X) # merging alpha-beta components for spin unrestricted representation #TODO: should be included into function not in main
     true_file = path+'/data/H2O_spahm_b.npy_alpha_beta.npy'
     X_true = np.load(true_file)
     assert(X_true.shape == X.shape)
     for Xa, Xa_true in zip(X, X_true):
         assert(np.linalg.norm(Xa-Xa_true) < 1e-8) # evaluating representation diff as norm (threshold = 1e-8)
 
+def test_water_closed():
+    path = os.path.dirname(os.path.realpath(__file__))
+    xyz_in = path+'/data/H2O.xyz'
+    mols = utils.load_mols([xyz_in], [None], [0], 'minao')
+    dms = utils.mols_guess(mols, [xyz_in], 'LB', spin=[None])
+    X = bond.get_repr(mols, [xyz_in], 'LB', spin=[None], with_symbols=False, same_basis=False)
+    true_file = path+'/data/H2O_spahm_b.npy'
+    X_true = np.load(true_file)
+    print(X_true.shape)
+    assert(X_true.shape == X.shape)
+    for Xa, Xa_true in zip(X, X_true):
+        assert(np.linalg.norm(Xa-Xa_true) < 1e-8) # evaluating representation diff as norm (threshold = 1e-8)
+
 def test_water_O_only():
     path = os.path.dirname(os.path.realpath(__file__))
     xyz_in = path+'/data/H2O.xyz'
     mols = utils.load_mols([xyz_in], [0], [0], 'minao')
     dms = utils.mols_guess(mols, [xyz_in], 'LB', spin=[0])
-    X = bond.bond(mols, dms, spin=[0], only_z=['O'])
+    X = bond.bond(mols, dms, only_z=['O'])
     X = np.squeeze(X) #contains a single elements but has shape (1,Nfeat)
     X = np.hstack(X) # merging alpha-beta components for spin unrestricted representation #TODO: should be included into function not in main
     true_file = path+'/data/H2O_spahm_b.npy_alpha_beta.npy'
@@ -36,7 +48,7 @@ def test_water_same_basis():
     xyz_in = path+'/data/H2O.xyz'
     mols = utils.load_mols([xyz_in], [0], [0], 'minao')
     dms = utils.mols_guess(mols, [xyz_in], 'LB', spin=[0])
-    X = bond.bond(mols, dms, spin=[0], same_basis=True)
+    X = bond.bond(mols, dms, same_basis=True)
     X = np.squeeze(X) #contains a single elements but has shape (1,Nfeat)
     X = np.hstack(X) # merging alpha-beta components for spin unrestricted representation #TODO: should be included into function not in main
     true_file = path+'/data/H2O_spahm_b_CCbas.npy_alpha_beta.npy'
@@ -48,9 +60,9 @@ def test_water_same_basis():
 def test_ecp():
     path = os.path.dirname(os.path.realpath(__file__))
     xyz_in = path+'/data/I2.xyz'
-    mols = utils.load_mols([xyz_in], [0], [None], 'minao', ecp='def2-svp')
-    dms = utils.mols_guess(mols, [xyz_in], 'LB', spin=[None])
-    X = bond.bond(mols, dms, spin=[None], same_basis=True)
+    mols = utils.load_mols([xyz_in], [0], [0], 'minao', ecp='def2-svp')
+    dms = utils.mols_guess(mols, [xyz_in], 'LB', spin=[0])
+    X = bond.bond(mols, dms, same_basis=True)
     X = np.squeeze(X) #contains a single elements but has shape (1,Nfeat)
     X = np.hstack(X) # merging alpha-beta components for spin unrestricted representation #TODO: should be included into function not in main
     true_file = path+'/data/I2_spahm-b_minao-def2-svp_alpha-beta.npy'
@@ -76,4 +88,5 @@ def test_from_list():
 if __name__ == '__main__':
     test_water()
     test_from_list()
+    test_water_closed()