small-fixes

environment.yml

@@ -3,51 +3,58 @@ channels:
   - conda-forge
   - defaults
 dependencies:
-  - _libgcc_mutex=0.1=main
-  - _openmp_mutex=4.5=1_gnu
-  - blas=1.0=mkl
-  - ca-certificates=2021.10.26=h06a4308_2
-  - certifi=2021.10.8=py39h06a4308_0
-  - intel-openmp=2021.4.0=h06a4308_3561
-  - ld_impl_linux-64=2.35.1=h7274673_9
-  - libffi=3.3=he6710b0_2
-  - libgcc-ng=9.3.0=h5101ec6_17
-  - libgomp=9.3.0=h5101ec6_17
-  - libstdcxx-ng=9.3.0=hd4cf53a_17
-  - mkl=2021.4.0=h06a4308_640
-  - mkl-service=2.4.0=py39h7f8727e_0
-  - mkl_fft=1.3.1=py39hd3c417c_0
-  - mkl_random=1.2.2=py39h51133e4_0
-  - ncurses=6.3=h7f8727e_2
-  - numpy=1.21.2=py39h20f2e39_0
-  - numpy-base=1.21.2=py39h79a1101_0
-  - openssl=1.1.1l=h7f8727e_0
-  - pip=22.3
-  - python=3.9.7=h12debd9_1
-  - readline=8.1=h27cfd23_0
-  - setuptools=58.0.4=py39h06a4308_0
-  - six=1.16.0=pyhd3eb1b0_0
-  - sqlite=3.37.0=hc218d9a_0
-  - tk=8.6.11=h1ccaba5_0
-  - tzdata=2021e=hda174b7_0
-  - wheel=0.37.0=pyhd3eb1b0_1
-  - xz=5.2.5=h7b6447c_0
-  - zlib=1.2.11=h7f8727e_4
-  - scikit-learn
+  - _libgcc_mutex=0.1=conda_forge
+  - _openmp_mutex=4.5=2_gnu
+  - bzip2=1.0.8=h4bc722e_7
+  - ca-certificates=2024.12.14=hbcca054_0
+  - ld_impl_linux-64=2.43=h712a8e2_2
+  - libexpat=2.6.4=h5888daf_0
+  - libffi=3.4.2=h7f98852_5
+  - libgcc=14.2.0=h77fa898_1
+  - libgcc-ng=14.2.0=h69a702a_1
+  - libgomp=14.2.0=h77fa898_1
+  - liblzma=5.6.3=hb9d3cd8_1
+  - libnsl=2.0.1=hd590300_0
+  - libsqlite=3.47.2=hee588c1_0
+  - libuuid=2.38.1=h0b41bf4_0
+  - libxcrypt=4.4.36=hd590300_1
+  - libzlib=1.3.1=hb9d3cd8_2
+  - ncurses=6.5=he02047a_1
+  - openssl=3.4.0=hb9d3cd8_0
+  - pip=24.3.1=pyh8b19718_2
+  - python=3.11.11=h9e4cc4f_1_cpython
+  - readline=8.2=h8228510_1
+  - setuptools=75.6.0=pyhff2d567_1
+  - tk=8.6.13=noxft_h4845f30_101
+  - tzdata=2024b=hc8b5060_0
+  - wheel=0.45.1=pyhd8ed1ab_1
   - pip:
-    - attrs==21.4.0
-    - h5py==3.6.0
-    - iniconfig==1.1.1
-    - packaging==21.3
-    - pluggy==1.0.0
-    - py==1.11.0
-    - pyparsing==3.0.6
-    - pyscf==2.0.1
-    - pytest==6.2.5
-    - scipy==1.7.3
-    - toml==0.10.2
-    - ase==3.22
-    - tqdm==4.66
-    - metatensor-core==0.1.8
-    - git+https://github.com/lcmd-epfl/cell2mol.git@22473bbf12a013467137a55a63c88fbbdc95baa2
-    - qstack/qstack-qml
+      - ase==3.22.0
+      - attrs==21.4.0
+      - git+https://github.com/lcmd-epfl/cell2mol.git
+      - certifi==2021.10.8
+      - contourpy==1.3.1
+      - cycler==0.12.1
+      - fonttools==4.55.3
+      - h5py==3.11.0
+      - iniconfig==1.1.1
+      - joblib==1.4.2
+      - kiwisolver==1.4.7
+      - matplotlib==3.10.0
+      - metatensor-core==0.1.8
+      - numpy==1.23.5
+      - packaging==21.3
+      - pillow==11.0.0
+      - pluggy==1.0.0
+      - py==1.11.0
+      - pyparsing==3.0.6
+      - pyscf==2.2.0
+      - pytest==6.2.5
+      - python-dateutil==2.9.0.post0
+      - qstack/qstack-qml
+      - scikit-learn==1.5.0
+      - scipy==1.10.0
+      - six==1.17.0
+      - threadpoolctl==3.5.0
+      - toml==0.10.2
+      - tqdm==4.66.0

qstack/spahm/rho/atom.py

@@ -63,7 +63,7 @@ def main():
     parser.add_argument('--aux-basis', dest='auxbasis',  default=defaults.auxbasis,            type=str, help=f"auxiliary basis set for density fitting (default: {defaults.auxbasis})")
     parser.add_argument('--model',     dest='model',     default=defaults.model,               type=str, help=f"the model to use when creating the representation (default: {defaults.model})")
     parser.add_argument('--dm',        dest='dm',        default=None,                         type=str, help="a density matrix to load instead of computing the guess")
-    parser.add_argument('--species',   dest='elements',  default=defaults.elements, nargs='+', type=str, help="the elements contained in the database")
+    parser.add_argument('--species',   dest='elements',  default=None, nargs='+', type=str, help="the elements contained in the database")
     parser.add_argument('--only',   dest='only_z', default=None, nargs='+', type=str, help="The restricted list of elements for which you want to generate the representation")
     parser.add_argument('--charge',    dest='charge',    default=0,                            type=int, help='total charge of the system (default: 0)')
     parser.add_argument('--spin',      dest='spin',      default=None,                         type=int, help='number of unpaired electrons (default: None) (use 0 to treat a closed-shell system in a UHF manner)')
@@ -77,7 +77,12 @@ def main():
     mol = compound.xyz_to_mol(check_file(args.mol), args.basis, charge=args.charge, spin=args.spin, unit=args.units, ecp=args.ecp)
     dm = None if args.dm is None else np.load(args.dm)
 
-    representations = get_repr(mol, args.elements, args.charge, args.spin,
+    if args.elements is None:
+        elements = sorted(mol.elements)
+    else:
+        elements = args.elements
+
+    representations = get_repr(mol, elements, args.charge, args.spin,
                                open_mod=args.omod,
                                dm=dm, guess=args.guess, model=args.model,
                                xc=args.xc, auxbasis=args.auxbasis, only_z=args.only_z)

qstack/spahm/rho/bond.py

@@ -3,13 +3,14 @@
 import os
 import argparse
 import numpy as np
+from itertools import chain
 from qstack.tools import correct_num_threads
 from . import utils, dmb_rep_bond as dmbb
 from .utils import defaults
 
 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.
 
@@ -19,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
@@ -39,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 spin is None:
-        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:
@@ -58,7 +56,7 @@ def bond(mols, dms,
     for imol, (mol, dm) in enumerate(zip(mols,dms)):
         if printlevel>0: print('mol', imol, flush=True)
         for iomod, omod in enumerate(omods):
-            DM  = utils.dm_open_mod(dm, omod) if spin else dm
+            DM  = utils.dm_open_mod(dm, omod)
             vec = dmbb.repr_for_mol(mol, DM, qqs, M, mybasis, idx, maxlen, cutoff, only_z=only_z)
             allvec[iomod,imol,:len(vec)] = vec
 
@@ -112,23 +110,25 @@ def get_repr(mols, xyzlist, guess,  xc=defaults.xc, spin=None, readdm=None,
         dms = []
 
     if len(only_z) > 0:
-        all_atoms   = np.array([z for mol in mols for z in mol.elements if z in only_z])
+        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]
-
     if split is False:
         shape  = (len(omods), -1, maxlen)
         atidx  = np.where(np.array([[1]*len(zin) + [0]*(natm-len(zin)) for zin in all_atoms]).flatten())
         allvec = allvec.reshape(shape)[:,atidx,:].reshape(shape)
-        all_atoms = all_atoms.flatten()
-        allvec = allvec.squeeze()
+        all_atoms = list(chain.from_iterable(all_atoms))
+        #allvec = allvec.squeeze()
     elif with_symbols:
         msg = f"You can not use 'split=True' and 'with_symbols=True' at the same time!"
         raise RuntimeError()
@@ -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)
@@ -192,10 +192,13 @@ def main():
     reps = get_repr(mols, xyzlist, args.guess, xc=args.xc, spin=spin, readdm=args.readdm, printlevel=args.print,
                       pairfile=args.pairfile, dump_and_exit=args.dump_and_exit, same_basis=args.same_basis,
                       bpath=args.bpath, cutoff=args.cutoff, omods=args.omod, with_symbols=args.with_symbols,
-                      elements=args.elements, only_m0=args.only_m0, zeros=args.zeros, split=args.split)
+                      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

@@ -21,21 +21,29 @@
 
 
 def get_chsp(fname, n):
-    if fname:
-        chsp = np.loadtxt(fname, dtype=int, ndmin=1)
+    def chsp_converter(chsp):
+        if chsp == 'None':
+            chsp = None
+        else:
+            chsp = int(chsp)
+        return chsp
+    if os.path.isfile(fname):
+        chsp = np.loadtxt(fname, dtype=object, converters=chsp_converter, encoding=None)
         if(len(chsp)!=n):
             raise RuntimeError(f'Wrong lengh of the file {fname}')
     else:
-        chsp = np.zeros(n, dtype=int)
+        raise RuntimeError(f"{fname} can not be found")
     return chsp
 
 
-def load_mols(xyzlist, charge, spin, basis, printlevel=0, units='ANG', ecp=None, progress=False):
+def load_mols(xyzlist, charge, spin, basis, printlevel=0, units='ANG', ecp=None, progress=False, srcdir=None):
     mols = []
     if progress:
         import tqdm
         xyzlist = tqdm.tqdm(xyzlist)
     for xyzfile, ch, sp in zip(xyzlist, charge, spin):
+        if srcdir is not None:
+            xyzfile = srcdir+xyzfile
         if printlevel>0: print(xyzfile, flush=True)
         mols.append(compound.xyz_to_mol(xyzfile, basis,
                                         charge=0 if ch is None else ch,
@@ -45,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:
@@ -66,7 +74,8 @@ def dm_open_mod(dm, omod):
     dmmod = {'sum':   lambda dm: dm[0]+dm[1],
              'diff':  lambda dm: dm[0]-dm[1],
              'alpha': lambda dm: dm[0],
-             'beta':  lambda dm: dm[1]}
+             'beta':  lambda dm: dm[1],
+             None:    lambda dm: dm}
     return dmmod[omod](dm)
 
 
@@ -134,8 +143,8 @@ def load_reps(f_in, from_list=True, srcdir=None, with_labels=False,
         else:
             is_single, is_labeled = file_format['is_single'], file_format['is_labeled']
         # if the given file contains a single representation create a one-element list
-        Xs = [np.load(f_in, allow_pickle=True)] if is_single else np.load(f_in, allow_pickle=True)
-    print(f"Loading {len(Xs)} representations (local = {local}, labeled = {is_labeled})")
+        Xs = [np.load(os.path.join(path2list,f_in), allow_pickle=True)] if is_single else np.load(os.path.join(path2list,f_in), allow_pickle=True)
+    if printlevel > 1: print(f"Loading {len(Xs)} representations (local = {local}, labeled = {is_labeled})")
     if progress:
         import tqdm
         Xs = tqdm.tqdm(Xs)
@@ -150,7 +159,10 @@ def load_reps(f_in, from_list=True, srcdir=None, with_labels=False,
                     reps.extend(x[:,1])
                     labels.extend(x[:,0])
             else:
-                reps.extend(x)
+                if sum_local:
+                    reps.append(x.sum(axis=0))
+                else:
+                    reps.extend(x)
         else:
            if is_labeled:
                 reps.append(x[1])

requirements.py3.11.txt

@@ -8,7 +8,7 @@ py==1.11.0
 pyparsing==3.0.6
 pyscf==2.2.0
 pytest==6.2.5
-numpy===1.22.3
+numpy==1.23.*
 scipy==1.10
 toml==0.10.2
 scikit-learn==1.5.0

tests/data/list_water.txt

@@ -0,0 +1,3 @@
+H2O_dist.xyz
+H2O.xyz
+rotated_H2O.xyz

tests/data/list_water_charges.txt

@@ -0,0 +1,3 @@
+0
+0
+0

tests/data/list_water_spins.txt

@@ -0,0 +1,3 @@
+0
+0
+0