Add LB for ECP

.gitignore

@@ -3,6 +3,9 @@ __pycache__/
 *.py[cod]
 *$py.class
 
+# vim tmp files
+*.swp
+
 # C extensions
 *.so
 

qstack/compound.py

@@ -10,15 +10,15 @@
 from qstack.tools import rotate_euler
 
 
-def xyz_to_mol(fin, basis="def2-svp", charge=0, spin=0, ignore=False, unit='ANG'):
+def xyz_to_mol(fin, basis="def2-svp", charge=0, spin=0, ignore=False, unit='ANG', ecp=None):
     """Reads a molecular file in xyz format and returns a pyscf Mole object.
-    
+
     Args:
         fin (str): Name (including path) of the xyz file to read.
         basis (str or dict): Basis set.
         charge (int): Charge of the molecule.
         spin (int): Spin of the molecule (alpha electrons - beta electrons).
-    
+
     Returns:
         A pyscf Mole object containing the molecule information.
     """
@@ -50,6 +50,9 @@ def xyz_to_mol(fin, basis="def2-svp", charge=0, spin=0, ignore=False, unit='ANG'
         mol.charge = -(sum(map(data.elements.charge, elements))%2)
         mol.spin = 0
 
+    if ecp:
+        mol.ecp = ecp
+
     mol.build()
     return mol
 
@@ -89,17 +92,17 @@ def mol_to_xyz(mol, fout, format="xyz"):
 
 
 def gmol_to_mol(fin, basis="def2-svp"):
-    """Reads a molecular file in gmol format and returns a pyscf Mole object.                                                     
-                                                                                                                                  
-    Args:                                                                                                                         
-        fin (str): Name (including path) of the xyz file to read.                                                                 
-        basis (str or dict): Basis set.                                                                                           
-        charge (int): Charge of the molecule.                                                                                     
-        spin (int): Alpha electrons - beta electrons.                                                                             
-                                                                                                                                  
-    Returns:                                                                                                                      
-        pyscf Mole: pyscf Mole object.                                                                                            
-    """             
+    """Reads a molecular file in gmol format and returns a pyscf Mole object.
+
+    Args:
+        fin (str): Name (including path) of the xyz file to read.
+        basis (str or dict): Basis set.
+        charge (int): Charge of the molecule.
+        spin (int): Alpha electrons - beta electrons.
+
+    Returns:
+        pyscf Mole: pyscf Mole object.
+    """
 
     from cell2mol.tmcharge_common import Cell, atom, molecule, ligand, metal
     from cell2mol.tmcharge_common import labels2formula
@@ -183,13 +186,13 @@ def gmol_to_mol(fin, basis="def2-svp"):
     return mole
 
 
-def make_auxmol(mol, basis):
+def make_auxmol(mol, basis, copy_ecp=False):
     """Builds an auxiliary Mole object given a basis set and a pyscf Mole object.
-    
+
     Args:
         mol (pyscf Mole): Original pyscf Mole object.
         basis (str or dict): Basis set.
-    
+
     Returns:
         An auxiliary pyscf Mole object.
     """
@@ -200,22 +203,24 @@ def make_auxmol(mol, basis):
     auxmol.charge = mol.charge
     auxmol.spin = mol.spin
     auxmol.basis = basis
+    if ecp:
+        auxmol.ecp = mol.ecp
     auxmol.build()
 
     return auxmol
 
 
 def rotate_molecule(mol, a, b, g, rad=False):
     """Rotate a molecule: transform nuclear coordinates given a set of Euler angles.
-    
+
     Args:
         mol (pyscf Mole): Original pyscf Mole object.
         a (float): Alpha Euler angle.
         b (float): Beta Euler angle.
         g (float): Gamma Euler angle.
         rad (bool) : Wheter the angles are in radians or not.
-    
-    
+
+
     Returns:
         A pyscf Mole object with transformed coordinates.
     """
@@ -237,10 +242,10 @@ def rotate_molecule(mol, a, b, g, rad=False):
 
 def fragments_read(frag_file):
     """Loads fragement definition from a frag file.
-    
+
     Args:
         frag_file (str): Name (including path) of the frag file to read.
-    
+
     Returns:
         A list of arrays containing the fragments.
     """
@@ -249,16 +254,16 @@ def fragments_read(frag_file):
     return fragments
 
 def fragment_partitioning(fragments, prop_atom_inp, normalize=True):
-    """Computes the contribution of each fragment.                                                                                
-                                                                                                                                  
-    Args:                                                                                                                         
-        fragments (numpy ndarray): Fragment definition                                                                            
-        prop_atom_inp (list of arrays or array): Coefficients densities.                                                          
-        normalize (bool): Normalized fragment partitioning. Defaults to True.                                                     
-                                                                                                                                  
-    Returns:                                                                                                                      
-        A list of arrays or an array containing the contribution of each fragment.                                                
-    """      
+    """Computes the contribution of each fragment.
+
+    Args:
+        fragments (numpy ndarray): Fragment definition
+        prop_atom_inp (list of arrays or array): Coefficients densities.
+        normalize (bool): Normalized fragment partitioning. Defaults to True.
+
+    Returns:
+        A list of arrays or an array containing the contribution of each fragment.
+    """
 
     if type(prop_atom_inp)==list:
         props_atom = prop_atom_inp

qstack/spahm/LB2020guess.py

@@ -1,6 +1,6 @@
 import copy
 import numpy
-import pyscf.data, pyscf.df
+import pyscf.data, pyscf.df, pyscf.scf
 
 """ Taken from https://github.com/briling/aepm and modified """
 
@@ -288,6 +288,30 @@ def use_charge(self, mol):
           acbasis = self.acbasis
       return acbasis
 
+  def use_ecp(self, mol, acbasis):
+      acbasis = copy.deepcopy(acbasis)
+      for iat, z in enumerate(mol.atom_charges()):
+          q = mol.atom_pure_symbol(iat)
+          zcore = pyscf.data.elements.charge(q) - z
+          if zcore > 0:
+              zrest = zcore
+              bad_idx = []
+              for iprim in range(len(acbasis[q])):
+                  if numpy.isclose(zrest, 0):
+                      break
+                  a, c = acbasis[q][iprim][1]
+                  renorm = self.renormalize(a)
+                  c /= renorm  # convert back to charge units: sum {c} == charge(q)
+                  dc = min(c, zrest)
+                  if numpy.isclose(c, dc):
+                      bad_idx.append(iprim)
+                  else:
+                      acbasis[q][iprim][1][1] = (c-dc)*renorm
+                  zrest -= dc
+              for i in bad_idx[::-1]:
+                  acbasis[q].pop(i)
+      return acbasis
+
   def get_auxweights(self, auxmol):
       w = numpy.zeros(auxmol.nao)
       iao = 0
@@ -307,17 +331,26 @@ def get_eri3c(self, mol, auxmol):
       eri3c      = pmol.intor('int3c2e_sph', shls_slice=shls_slice)
       return eri3c
 
+  def check_coefficients(self, mol, acbasis):
+      ch1 = sum(sum(c/self.renormalize(a) for _, (a, c) in acbasis[mol.atom_pure_symbol(iat)]) for iat in range(mol.natm))
+      ch2 = sum(mol.atom_charges()) - (mol.charge if self.parameters == 'HF' else 0)
+      assert numpy.isclose(ch1, ch2)
+
   def HLB20(self, mol):
       acbasis = self.use_charge(mol)
+      if mol.has_ecp():
+          acbasis = self.use_ecp(mol, acbasis)
+      self.check_coefficients(mol, acbasis)
+
       auxmol  = pyscf.df.make_auxmol(mol, acbasis)
       eri3c   = self.get_eri3c(mol, auxmol)
       auxw    = self.get_auxweights(auxmol)
       return self.merge_caps(auxw, eri3c)
 
   def Heff(self, mol):
-      self.mol = mol
-      self.Hcore = mol.intor('int1e_nuc_sph') + mol.intor('int1e_kin_sph')
-      self.H    = self.Hcore + self.HLB20(mol)
+      self.mol   = mol
+      self.Hcore = pyscf.scf.hf.get_hcore(mol)
+      self.H     = self.Hcore + self.HLB20(mol)
       return self.H
 
 
@@ -340,6 +373,9 @@ def HLB20_ints_deriv(iat):
 
   def HLB20_generator(self, mol):
       acbasis = self.use_charge(mol)
+      if mol.has_ecp():
+          acbasis = self.use_ecp(mol, acbasis)
+      self.check_coefficients(mol, acbasis)
       auxmol  = pyscf.df.make_auxmol(mol, acbasis)
       eri3c   = self.HLB20_ints_generator(mol, auxmol)
       auxw    = self.get_auxweights(auxmol)

qstack/spahm/compute_spahm.py

@@ -3,12 +3,12 @@
 
 def get_guess_orbitals(mol, guess, xc="pbe"):
     """ Compute the guess Hamiltonian.
-    
+
     Args:
         mol (pyscf Mole): pyscf Mole object.
         guess (funct): Method used to compute the guess Hamiltonian. Output of get_guess.
         xc (str): Exchange-correlation functional. Defaults to pbe.
-    
+
     Returns:
         A 1D numpy array containing the eigenvalues and a 2D numpy array containing the eigenvectors of the guess Hamiltonian.
     """
@@ -28,13 +28,13 @@ def get_guess_orbitals_grad(mol, guess):
 
 def get_guess_dm(mol, guess, xc="pbe", openshell=None):
     """ Compute the density matrix with the guess Hamiltonian.
-    
+
     Args:
         mol (pyscf Mole): pyscf Mole object.
         guess (funct): Method used to compute the guess Hamiltonian. Output of get_guess.
         xc (str): Exchange-correlation functional. Defaults to pbe
         openshell (bool): . Defaults to None.
-    
+
     Returns:
         A numpy ndarray containing the density matrix computed using the guess Hamiltonian.
     """
@@ -43,12 +43,12 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None):
 
 def get_spahm_representation(mol, guess_in, xc="pbe"):
     """ Compute the SPAHM representation.
-    
+
     Args:
         mol (pyscf Mole): pyscf Mole object.
         guess_in (str): Method used to obtain the guess Hamiltoninan.
         xc (str): Exchange-correlation functional. Defaults to pbe.
-    
+
     Returns:
         A numpy ndarray containing the SPAHM representation.
     """

qstack/spahm/rho/utils.py

@@ -29,11 +29,14 @@ def get_chsp(fname, n):
     return chsp
 
 
-def load_mols(xyzlist, charge, spin, basis, printlevel=0, units='ANG'):
+def load_mols(xyzlist, charge, spin, basis, printlevel=0, units='ANG', ecp=None):
     mols = []
     for xyzfile, ch, sp in zip(xyzlist, charge, spin):
         if printlevel>0: print(xyzfile, flush=True)
-        mols.append(compound.xyz_to_mol(xyzfile, basis, charge=0 if ch is None else ch, spin=0 if ch is None else sp, unit=units)) #TODO
+        mols.append(compound.xyz_to_mol(xyzfile, basis,
+                                        charge=0 if ch is None else ch,
+                                        spin=0 if ch is None else sp,
+                                        unit=units, ecp=ecp)) #TODO
     if printlevel>0: print()
     return mols
 

tests/data/H2Te.xyz

@@ -0,0 +1,5 @@
+3
+
+ Te     0.000000    0.000000   -0.795911
+ H     -1.186112    0.000000    0.397955
+ H      1.186112    0.000000    0.397955

tests/test_spahm.py

@@ -5,6 +5,7 @@
 from qstack import compound
 from qstack.spahm import compute_spahm
 
+
 def test_spahm_huckel():
     path = os.path.dirname(os.path.realpath(__file__))
     mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'def2svp', charge=0, spin=0)
@@ -14,6 +15,7 @@ def test_spahm_huckel():
     assert(R.shape == (2,5))
     assert(np.abs(np.sum(R-true_R)) < 1e-05)
 
+
 def test_spahm_LB():
     path = os.path.dirname(os.path.realpath(__file__))
     mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'def2svp', charge=1, spin=1)
@@ -23,6 +25,19 @@ def test_spahm_LB():
     assert(R.shape == (2,5))
     assert(np.abs(np.sum(R-true_R)) < 1e-05)
 
+
+def test_spahm_LB_ecp():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2Te.xyz', 'minao',
+                              ecp='def2-svp', charge=0, spin=0)
+    R = compute_spahm.get_spahm_representation(mol, 'lb')[0]
+    true_R = np.array([-5.68297474, -3.91180966, -3.91176332, -3.90721427, -1.22967252, -1.22469672,
+                       -1.22145412, -1.2210437 , -1.22099792, -0.43285024, -0.20943343, -0.15915716,
+                       -0.07260264])
+    assert np.allclose(R, true_R)
+
+
 if __name__ == '__main__':
     test_spahm_huckel()
     test_spahm_LB()
+    test_spahm_LB_ecp()

tests/test_spahm_grad.py

@@ -7,16 +7,19 @@
 
 np.set_printoptions(formatter={'all':lambda x: '% .4f'%x})
 
+
 def build_mol(mol, r):
     r = r.reshape((-1,3))
     mymol = gto.Mole()
     mymol.charge = mol.charge
     mymol.spin   = mol.spin
+    mymol.ecp    = mol.ecp
     mymol.atom = [ (mol.atom_symbol(i), r[i]) for i in range(mol.natm)]
     mymol.basis = mol.basis
     mymol.build()
     return mymol
 
+
 def grad_num(func, mol, guess, eps=1e-4):
     r = mol.atom_coords(unit='ang').flatten()
     g = []
@@ -29,6 +32,7 @@ def grad_num(func, mol, guess, eps=1e-4):
         g.append((8.0*e1-8.0*e2 + e22-e11) / (12.0*eps))
     return np.array(g)*data.nist.BOHR
 
+
 def test_spahm_ev_grad():
     def spahm_ev(r, mol, guess):
         mymol = build_mol(mol, r)
@@ -42,6 +46,7 @@ def spahm_ev(r, mol, guess):
     for g1, g2 in zip(ngrad.T, agrad.reshape(-1, 9)):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
+
 def test_spahm_re_grad():
     def spahm_re(r, mol, guess_in):
         mymol = build_mol(mol, r)
@@ -56,6 +61,21 @@ def spahm_re(r, mol, guess_in):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
+def test_spahm_ev_grad_ecp():
+    def spahm_ev(r, mol, guess):
+        mymol = build_mol(mol, r)
+        e, c = spahm.compute_spahm.get_guess_orbitals(mymol, guess[0])
+        return e
+    path  = os.path.dirname(os.path.realpath(__file__))
+    mol   = compound.xyz_to_mol(path+'/data/H2Te.xyz', 'minao', charge=0, spin=0, ecp='def2-svp')
+    guess = spahm.guesses.get_guess_g('lb')
+    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess)
+    ngrad = grad_num(spahm_ev, mol, guess)
+    for g1, g2 in zip(ngrad.T, agrad.reshape(-1, 9)):
+        assert(np.linalg.norm(g1-g2)<1e-6)
+
+
 if __name__ == '__main__':
     test_spahm_ev_grad()
     test_spahm_re_grad()
+    test_spahm_ev_grad_ecp()