Merge branch 'master' into mcfun-pbc

CHANGELOG

@@ -1,3 +1,8 @@
+PySCF (2023-08-01)
+-----------------------
+* Added
+  - NVT Molecular Dynamics
+
 PySCF 2.3.0 (2023-07-04)
 ------------------------
 * Added

FEATURES

@@ -206,7 +206,7 @@
     Libcint library
   - Support for 4-index integral transformation with 4 different orbitals
   - PBC 2-electron MO integrals
-  - Integral transformation for (4-component and 2-compoent relativistic) spinor
+  - Integral transformation for (4-component and 2-component relativistic) spinor
     integrals
 
 

examples/1-advanced/033-constrained_dft.py

@@ -328,7 +328,7 @@ def get_fock(h1e, s1e, vhf, dm, cycle=0, mf_diis=None):
             elif diis_type == 3:
                 fock = cdft_diis.update(fock, scf.diis.get_err_vec(s1e, dm, fock))
             else:
-                print("\nWARN: Unknow CDFT DIIS type, NO DIIS IS USED!!!\n")
+                print("\nWARN: Unknown CDFT DIIS type, NO DIIS IS USED!!!\n")
 
         if diis_pos == 'post' or diis_pos == 'both':
             cdft_conv_flag = False

examples/dft/11-grid_scheme.py

@@ -49,7 +49,7 @@
 # Stratmann-Scuseria weight scheme
 method = dft.RKS(mol)
 method.grids.becke_scheme = dft.stratmann
-print('Changed grid partition funciton.  E = %.12f' % method.kernel())
+print('Changed grid partition function.  E = %.12f' % method.kernel())
 
 # Grids level 0 - 9.  Big number indicates dense grids. Default is 3
 method = dft.RKS(mol)
@@ -71,4 +71,4 @@
 #grids.prune = dft.sg1_prune
 method = dft.RKS(mol)
 method.grids.prune = None
-print('Changed grid partition funciton.  E = %.12f' % method.kernel())
+print('Changed grid partition function.  E = %.12f' % method.kernel())

examples/local_orb/pmloc.py

@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 #
-# Localization funciton: pmloc.loc(mol,mocoeff)
+# Localization function: pmloc.loc(mol,mocoeff)
 #
 # -ZL@20151227- Add Boys. Note that to view LMOs
 #		in jmol, 'orbital energies' need

examples/md/00-simple_nve.py

@@ -24,13 +24,13 @@
 myintegrator = pyscf.md.NVE(myscanner,
                             dt=5,
                             steps=10,
-                            energy_output="BOMD.md.energies",
+                            data_output="BOMD.md.data",
                             trajectory_output="BOMD.md.xyz").run()
 
 # Note that we can also just pass the CASSCF object directly to
 # generate the integrator and it will automatically convert it to a scanner
 # myintegrator = pyscf.md.NVE(mycas, dt=5, steps=100)
 
 # Close the file streams for the energy and trajectory.
-myintegrator.energy_output.close()
+myintegrator.data_output.close()
 myintegrator.trajectory_output.close()

examples/md/03-mb_initial_veloc.py

@@ -37,7 +37,7 @@
 myhf = mol.RHF()
 
 # We set the initial velocity by passing to "veloc"
-myintegrator = pyscf.md.NVE(myhf, dt=5, steps=10, veloc=init_veloc)
+myintegrator = pyscf.md.NVE(myhf, dt=5, steps=10, veloc=init_veloc, data_output="NVE.md.data")
 
 myintegrator.run()
 

examples/md/04-mb_nvt_berendson.py

@@ -0,0 +1,29 @@
+#!/usr/bin/env python
+
+'''
+Run an NVT BOMD simulation using initial velocity from 
+a Maxwell-Boltzmann distribution at 300K. 
+'''
+
+import pyscf
+import pyscf.md
+
+mol = pyscf.M(
+    atom='''O 0 0 0; H 0 -0.757 0.587; H 0 0.757 0.587''',
+    basis = 'ccpvdz')
+
+myhf = mol.RHF()
+
+# initial velocities from a Maxwell-Boltzmann distribution [T in K and velocities are returned in (Bohr/ time a.u.)]
+init_veloc = pyscf.md.distributions.MaxwellBoltzmannVelocity(mol, T=300)
+
+# We set the initial velocity by passing to "veloc", 
+#T is the ensemble temperature in K and taut is the Berendsen Thermostat time constant given in time a.u.
+myintegrator = pyscf.md.integrators.NVTBerendson(myhf, dt=5, steps=100, 
+			     T=300, taut=50, veloc=init_veloc,
+			     data_output="NVT.md.data", 
+			     trajectory_output="NVT.md.xyz").run()
+
+# Close the file streams for the energy, temperature and trajectory.
+myintegrator.data_output.close()
+myintegrator.trajectory_output.close()

examples/scf/32-break_spin_symm.py

@@ -29,7 +29,7 @@
 # (alpha,beta).  Assigning alpha-DM and beta-DM to different value can break
 # the spin symmetry.
 #
-# In the following example, the funciton get_init_guess returns the
+# In the following example, the function get_init_guess returns the
 # superposition of atomic density matrices in which the alpha and beta
 # components are degenerated.  The degeneracy are destroyed by zeroing out the
 # beta 1s,2s components.

examples/symm/33-lz_adaption.py

@@ -40,7 +40,7 @@ def lz_adapt_(mol):
         mol.symm_orb[Ex] = lz_plus   # x+iy
     return mol
 
-mol = gto.M(atom='Ne', basis='ccpvtz', symmetry=True)
+mol = gto.M(atom='Ne', basis='ccpvtz', symmetry='d2h')
 mol = lz_adapt_(mol)
 mf = scf.RHF(mol)
 mf.run()

pyscf/df/addons.py

@@ -29,6 +29,7 @@
 ETB_BETA = getattr(__config__, 'df_addons_aug_dfbasis', 2.0)
 FIRST_ETB_ELEMENT = getattr(__config__, 'df_addons_aug_start_at', 36)  # 'Rb'
 
+
 # Obtained from http://www.psicode.org/psi4manual/master/basissets_byfamily.html
 DEFAULT_AUXBASIS = {
     # AO basis       JK-fit                     MP2-fit

pyscf/df/incore.py

@@ -31,7 +31,7 @@
 LINEAR_DEP_THR = getattr(__config__, 'df_df_DF_lindep', 1e-12)
 
 
-# This funciton is aliased for backward compatibility.
+# This function is aliased for backward compatibility.
 format_aux_basis = addons.make_auxmol
 
 

pyscf/dft/libxc.py

@@ -708,12 +708,12 @@ def available_libxc_functionals():
 def _xc_key_without_underscore(xc_keys):
     new_xc = []
     for key, xc_id in xc_keys.items():
-        for delimeter in ('_XC_', '_X_', '_C_', '_K_'):
-            if delimeter in key:
-                key0, key1 = key.split(delimeter)
+        for delimiter in ('_XC_', '_X_', '_C_', '_K_'):
+            if delimiter in key:
+                key0, key1 = key.split(delimiter)
                 new_key1 = key1.replace('_', '').replace('-', '')
                 if key1 != new_key1:
-                    new_xc.append((key0+delimeter+new_key1, xc_id))
+                    new_xc.append((key0+delimiter+new_key1, xc_id))
                 break
     return new_xc
 XC_CODES.update(_xc_key_without_underscore(XC_CODES))
@@ -1060,7 +1060,7 @@ def parse_xc(description):
         part blank. E.g. description='slater,' means pure LDA functional.
       - To neglect X functional (just apply C functional), leave the first
         part blank. E.g. description=',vwn' means pure VWN functional.
-      - If compound XC functional is specified, no matter whehter it is in the
+      - If compound XC functional is specified, no matter whether it is in the
         X part (the string in front of comma) or the C part (the string behind
         comma), both X and C functionals of the compound XC functional will be
         used.
@@ -1485,7 +1485,7 @@ def _eval_xc(hyb, fn_facs, rho, spin=0, relativity=0, deriv=1, verbose=None):
         for rho_ud in [rho_u, rho_d]:
             assert rho_ud.shape[0] >= 6
     else:
-        raise ValueError("Unknow nvar {}".format(nvar))
+        raise ValueError("Unknown nvar {}".format(nvar))
 
     outlen = (math.factorial(nvar+deriv) //
               (math.factorial(nvar) * math.factorial(deriv)))

pyscf/dft/numint.py

@@ -86,7 +86,7 @@ def eval_ao(mol, coords, deriv=0, shls_slice=None,
         2D array of shape (N,nao) for AO values if deriv = 0.
         Or 3D array of shape (:,N,nao) for AO values and AO derivatives if deriv > 0.
         In the 3D array, the first (N,nao) elements are the AO values,
-        followed by (3,N,nao) for x,y,z compoents;
+        followed by (3,N,nao) for x,y,z components;
         Then 2nd derivatives (6,N,nao) for xx, xy, xz, yy, yz, zz;
         Then 3rd derivatives (10,N,nao) for xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz;
         ...
@@ -351,7 +351,7 @@ def eval_rho2(mol, ao, mo_coeff, mo_occ, non0tab=None, xctype='LDA',
         xctype : str
             LDA/GGA/mGGA.  It affects the shape of the return density.
         with_lapl: bool
-            Wether to compute laplacian. It affects the shape of returns.
+            Whether to compute laplacian. It affects the shape of returns.
         verbose : int or object of :class:`Logger`
             No effects.
 

pyscf/dft/rks.py

@@ -199,7 +199,7 @@ def _get_k_lr(mol, dm, omega=0, hermi=0, vhfopt=None):
                      'It is replaced by mol.get_k(mol, dm, omege=omega)')
     dm = numpy.asarray(dm)
 # Note, ks object caches the ERIs for small systems. The cached eris are
-# computed with regular Coulomb operator. ks.get_jk or ks.get_k do not evalute
+# computed with regular Coulomb operator. ks.get_jk or ks.get_k do not evaluate
 # the K matrix with the range separated Coulomb operator.  Here jk.get_jk
 # function computes the K matrix with the modified Coulomb operator.
     nao = dm.shape[-1]

pyscf/eph/rhf.py

@@ -39,7 +39,7 @@ def kernel(ephobj, mo_energy=None, mo_coeff=None, mo_occ=None, mo_rep=False):
 
     # chkfile is used to pass first orbitals from hessian methods to eph methods
     # TODO: Remove the dependence to chfile and return first orbitals from a function
-    assert ephobj.chkfile is not None, 'chkfile is requred to save first order orbitals'
+    assert ephobj.chkfile is not None, 'chkfile is required to save first order orbitals'
 
     de = ephobj.hess_elec(mo_energy, mo_coeff, mo_occ)
     ephobj.de = de + ephobj.hess_nuc(ephobj.mol)

pyscf/fci/cistring.py

@@ -32,7 +32,7 @@ def make_strings(orb_list, nelec):
         orbital, bit-1 means occupied and bit-0 means unoccupied.  The lowest
         (right-most) bit corresponds to the lowest orbital in the orb_list.
 
-    Exampels:
+    Examples:
 
     >>> [bin(x) for x in make_strings((0,1,2,3),2)]
     [0b11, 0b101, 0b110, 0b1001, 0b1010, 0b1100]

pyscf/fci/direct_spin1_cyl_sym.py

@@ -26,7 +26,7 @@
 Hamiltonian. For 2D irreps, the real part and the imaginary part of the complex
 FCI wavefunction are identical to the Ex and Ey wavefunction obtained from
 direct_spin1_symm module. However, any observables from the complex FCI
-wavefunction should have an indentical one from either Ex or Ey wavefunction
+wavefunction should have an identical one from either Ex or Ey wavefunction
 of direct_spin1_symm.
 '''
 

pyscf/fci/direct_spin1_symm.py

@@ -414,7 +414,7 @@ def reorder_eri(eri, norb, orbsym):
 
     # irrep of (ij| pair
     trilirrep = (orbsym[:,None] ^ orbsym)[np.tril_indices(norb)]
-    # and the number of occurence for each irrep
+    # and the number of occurrences for each irrep
     dimirrep = np.asarray(np.bincount(trilirrep), dtype=np.int32)
     # we sort the irreps of (ij| pair, to group the pairs which have same irreps
     # "order" is irrep-id-sorted index. The (ij| paired is ordered that the

pyscf/fci/rdm.py

@@ -101,8 +101,8 @@ def make_rdm1_spin1(fname, cibra, ciket, norb, nelec, link_index=None):
         link_indexa, link_indexb = link_index
     na,nlinka = link_indexa.shape[:2]
     nb,nlinkb = link_indexb.shape[:2]
-    assert (cibra.size == na*nb)
-    assert (ciket.size == na*nb)
+    assert (cibra.size == na*nb), '{} {} {}'.format (cibra.size, na, nb)
+    assert (ciket.size == na*nb), '{} {} {}'.format (ciket.size, na, nb)
     rdm1 = numpy.empty((norb,norb))
     fn = getattr(librdm, fname)
     fn(rdm1.ctypes.data_as(ctypes.c_void_p),

pyscf/fci/selected_ci.py

@@ -902,7 +902,7 @@ def _all_linkstr_index(ci_strs, norb, nelec):
     dd_indexb = des_des_linkstr_tril(ci_strs[1], norb, nelec[1])
     return cd_indexa, dd_indexa, cd_indexb, dd_indexb
 
-# numpy.ndarray does not allow to attach attribtues.  Overwrite the
+# numpy.ndarray does not allow to attach attributes.  Overwrite the
 # numpy.ndarray class to tag the ._strs attribute
 class SCIvector(numpy.ndarray):
     '''An 2D np array for selected CI coefficients'''

pyscf/fci/selected_ci_symm.py

@@ -38,7 +38,7 @@ def reorder4irrep(eri, norb, link_index, orbsym, offdiag=0):
     orbsym = orbsym % 10
     # irrep of (ij| pair
     trilirrep = (orbsym[:,None] ^ orbsym)[numpy.tril_indices(norb, offdiag)]
-    # and the number of occurence for each irrep
+    # and the number of occurrences for each irrep
     dimirrep = numpy.array(numpy.bincount(trilirrep), dtype=numpy.int32)
     # we sort the irreps of (ij| pair, to group the pairs which have same irreps
     # "order" is irrep-id-sorted index. The (ij| paired is ordered that the

pyscf/fci/spin_op.py

@@ -68,7 +68,7 @@ def spin_square_general(dm1a, dm1b, dm2aa, dm2ab, dm2bb, mo_coeff, ovlp=1):
                          + Gamma_{i\beta k\beta ,j\beta l\beta})
                          + (n_\alpha+n_\beta)/4
 
-    Given the overlap betwen non-degenerate alpha and beta orbitals, this
+    Given the overlap between non-degenerate alpha and beta orbitals, this
     function can compute the expectation value spin square operator for
     UHF-FCI wavefunction
     '''

pyscf/geomopt/berny_solver.py

@@ -131,12 +131,12 @@ def kernel(method, assert_convergence=ASSERT_CONV,
     # When symmetry is enabled, the molecule may be shifted or rotated to make
     # the z-axis be the main axis. The transformation can cause inconsistency
     # between the optimization steps. The transformation is muted by setting
-    # an explict point group to the keyword mol.symmetry (see symmetry
+    # an explicit point group to the keyword mol.symmetry (see symmetry
     # detection code in Mole.build function).
     if mol.symmetry:
         mol.symmetry = mol.topgroup
 
-# temporary interface, taken from berny.py optimize function
+    # temporary interface, taken from berny.py optimize function
     berny_log = to_berny_log(log)
     geom = to_berny_geom(mol, include_ghost)
     optimizer = Berny(geom, logger=berny_log, **kwargs)

pyscf/geomopt/geometric_solver.py

@@ -132,13 +132,13 @@ def kernel(method, assert_convergence=ASSERT_CONV,
     engine = PySCFEngine(g_scanner)
     engine.callback = callback
     engine.maxsteps = maxsteps
-    # To avoid overwritting method.mol
+    # To avoid overwriting method.mol
     engine.mol = g_scanner.mol.copy()
 
     # When symmetry is enabled, the molecule may be shifted or rotated to make
     # the z-axis be the main axis. The transformation can cause inconsistency
     # between the optimization steps. The transformation is muted by setting
-    # an explict point group to the keyword mol.symmetry (see symmetry
+    # an explicit point group to the keyword mol.symmetry (see symmetry
     # detection code in Mole.build function).
     if engine.mol.symmetry:
         engine.mol.symmetry = engine.mol.topgroup

pyscf/grad/lagrange.py

@@ -118,7 +118,7 @@ def my_Lvec_last ():
             return Lvec_last
         precond = self.get_lagrange_precond (Adiag, level_shift=level_shift, **kwargs)
         it = np.asarray ([0])
-        logger.debug(self, 'Lagrange multiplier determination intial gradient norm: %.8g',
+        logger.debug(self, 'Lagrange multiplier determination initial gradient norm: %.8g',
                      linalg.norm(bvec))
         my_call = self.get_lagrange_callback (Lvec_last, it, my_geff)
         Aop_obj = sparse_linalg.LinearOperator ((self.nlag,self.nlag), matvec=Aop,

pyscf/grad/rhf.py

@@ -93,21 +93,19 @@ def grad_nuc(mol, atmlst=None):
     '''
     Derivatives of nuclear repulsion energy wrt nuclear coordinates
     '''
-    gs = numpy.zeros((mol.natm,3))
-    for j in range(mol.natm):
-        q2 = mol.atom_charge(j)
-        r2 = mol.atom_coord(j)
-        for i in range(mol.natm):
-            if i != j:
-                q1 = mol.atom_charge(i)
-                r1 = mol.atom_coord(i)
-                r = numpy.sqrt(numpy.dot(r1-r2,r1-r2))
-                gs[j] -= q1 * q2 * (r2-r1) / r**3
+    z = mol.atom_charges()
+    r = mol.atom_coords()
+    dr = r[:,None,:] - r
+    dist = numpy.linalg.norm(dr, axis=2)
+    diag_idx = numpy.diag_indices(z.size)
+    dist[diag_idx] = 1e100
+    rinv = 1./dist
+    rinv[diag_idx] = 0.
+    gs = numpy.einsum('i,j,ijx,ij->ix', -z, z, dr, rinv**3)
     if atmlst is not None:
         gs = gs[atmlst]
     return gs
 
-
 def get_hcore(mol):
     '''Part of the nuclear gradients of core Hamiltonian'''
     h = mol.intor('int1e_ipkin', comp=3)

pyscf/grad/test/test_rhf.py

@@ -205,8 +205,25 @@ def test_grad_with_symmetry(self):
         g_x = scf.RHF (mol).run ().nuc_grad_method ().kernel ()
         self.assertAlmostEqual(abs(ref[:,2] - g_x[:,0]).max(), 0, 9)
 
+    def test_grad_nuc(self):
+        mol = gto.M(atom='He 0 0 0; He 0 1 2; H 1 2 1; H 1 0 0')
+        gs = grad.rhf.grad_nuc(mol)
+        ref = grad_nuc(mol)
+        self.assertAlmostEqual(abs(gs - ref).max(), 0, 9)
+
+def grad_nuc(mol):
+    gs = numpy.zeros((mol.natm,3))
+    for j in range(mol.natm):
+        q2 = mol.atom_charge(j)
+        r2 = mol.atom_coord(j)
+        for i in range(mol.natm):
+            if i != j:
+                q1 = mol.atom_charge(i)
+                r1 = mol.atom_coord(i)
+                r = numpy.sqrt(numpy.dot(r1-r2,r1-r2))
+                gs[j] -= q1 * q2 * (r2-r1) / r**3
+    return gs
 
 if __name__ == "__main__":
     print("Full Tests for RHF Gradients")
     unittest.main()
-

pyscf/gw/rpa.py

@@ -74,7 +74,7 @@ def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=40, x0=0.5, verbose=logger.NOT
     # Compute RPA correlation energy
     e_corr = get_rpa_ecorr(rpa, Lpq, freqs, wts)
 
-    # Compute totol energy
+    # Compute total energy
     e_tot = e_hf + e_corr
 
     logger.debug(rpa, '  RPA total energy = %s', e_tot)
@@ -232,7 +232,7 @@ def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40, x0=0.5):
             mo_energy : 1D array (nmo), mean-field mo energy
             mo_coeff : 2D array (nmo, nmo), mean-field mo coefficient
             Lpq : 3D array (naux, nmo, nmo), 3-index ERI
-            nw: interger, grid number
+            nw: integer, grid number
             x0: real, scaling factor for frequency grid
 
         Returns: