From 0c6a0aa3ea2d0cc8193c9d708502abcc477cd4c9 Mon Sep 17 00:00:00 2001 From: Aniruddha Seal Date: Fri, 15 Sep 2023 13:00:16 -0500 Subject: [PATCH 1/6] added NVT BOMD and example (#1818) * add NVT BOMD and example * overload functions * remove default NVT * remove NVTBerendson * update init * update example * add NVTBerendson unit test * corrected Temperature() * corrected temperature * added temp rmsf check * added temperature test * removed pandas dependency * added temperature test linear mol * corrected flake8 errors --- CHANGELOG | 5 + examples/md/00-simple_nve.py | 4 +- examples/md/03-mb_initial_veloc.py | 2 +- examples/md/04-mb_nvt_berendson.py | 29 +++++ pyscf/md/__init__.py | 1 - pyscf/md/integrators.py | 166 ++++++++++++++++++++++++----- pyscf/md/test/test_NVTBerendson.py | 72 +++++++++++++ pyscf/md/test/test_prop.py | 89 ++++++++++++++++ 8 files changed, 336 insertions(+), 32 deletions(-) create mode 100644 examples/md/04-mb_nvt_berendson.py create mode 100644 pyscf/md/test/test_NVTBerendson.py create mode 100644 pyscf/md/test/test_prop.py diff --git a/CHANGELOG b/CHANGELOG index 99d40ac9b1..3357410bef 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,3 +1,8 @@ +PySCF (2023-08-01) +----------------------- +* Added + - NVT Molecular Dynamics + PySCF 2.3.0 (2023-07-04) ------------------------ * Added diff --git a/examples/md/00-simple_nve.py b/examples/md/00-simple_nve.py index 611f1ae6ed..0f5ae330b6 100644 --- a/examples/md/00-simple_nve.py +++ b/examples/md/00-simple_nve.py @@ -24,7 +24,7 @@ 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 @@ -32,5 +32,5 @@ # 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() diff --git a/examples/md/03-mb_initial_veloc.py b/examples/md/03-mb_initial_veloc.py index 46f79bd121..c209d2483b 100644 --- a/examples/md/03-mb_initial_veloc.py +++ b/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() diff --git a/examples/md/04-mb_nvt_berendson.py b/examples/md/04-mb_nvt_berendson.py new file mode 100644 index 0000000000..7a0a7dabdd --- /dev/null +++ b/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() diff --git a/pyscf/md/__init__.py b/pyscf/md/__init__.py index 688ffb0f27..72ffe99d58 100644 --- a/pyscf/md/__init__.py +++ b/pyscf/md/__init__.py @@ -45,4 +45,3 @@ def set_seed(seed): from pyscf.md import integrators, distributions NVE = integrators.VelocityVerlet - diff --git a/pyscf/md/integrators.py b/pyscf/md/integrators.py index e9f5ceabf2..a303331e1c 100644 --- a/pyscf/md/integrators.py +++ b/pyscf/md/integrators.py @@ -13,7 +13,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # -# Author: Matthew Hennefarth +# Authors: Matthew Hennefarth , +# Aniruddha Seal import os import numpy as np @@ -156,8 +157,9 @@ class _Integrator(lib.StreamObject): stdout : file object Default is self.scanner.mol.stdout. - energy_output : file object - Stream to write energy to during the course of the simulation. + data_output : file object + Stream to write energy and temperature to + during the course of the simulation. trajectory_output : file object Stream to write the trajectory to during the course of the @@ -206,7 +208,7 @@ def __init__(self, method, **kwargs): self.epot = None self.ekin = None self.time = 0 - self.energy_output = None + self.data_output = None self.trajectory_output = None self.callback = None @@ -239,7 +241,8 @@ def kernel(self, veloc=None, steps=None, dump_flags=True, verbose=None): Print level Returns: - _Integrator with final epot, ekin, mol, and veloc of the simulation. + _Integrator with final epot, ekin, temp, + mol, and veloc of the simulation. ''' if veloc is not None: @@ -261,23 +264,23 @@ def kernel(self, veloc=None, steps=None, dump_flags=True, verbose=None): data.elements.COMMON_ISOTOPE_MASSES[m] * data.nist.AMU2AU for m in self.mol.atom_charges()]) - # avoid opening energy_output file twice - if type(self.energy_output) is str: + # avoid opening data_output file twice + if type(self.data_output) is str: if self.verbose > logger.QUIET: - if os.path.isfile(self.energy_output): - print('overwrite energy output file: %s' % - self.energy_output) + if os.path.isfile(self.data_output): + print('overwrite data output file: %s' % + self.data_output) else: - print('energy output file: %s' % self.energy_output) + print('data output file: %s' % self.data_output) - if self.energy_output == '/dev/null': - self.energy_output = open(os.devnull, 'w') + if self.data_output == '/dev/null': + self.data_output = open(os.devnull, 'w') else: - self.energy_output = open(self.energy_output, 'w') - self.energy_output.write( + self.data_output = open(self.data_output, 'w') + self.data_output.write( 'time Epot Ekin ' - 'Etot\n' + 'Etot T\n' ) # avoid opening trajectory_output file twice @@ -337,12 +340,13 @@ def compute_kinetic_energy(self): def temperature(self): '''Returns the temperature of the system''' + # checked against ORCA for linear and non-linear molecules dof = 3 * len(self.mol.atom_coords()) # Temp = 2/(3*k*N_f) * KE # = 2/(3*k*N_f)*\sum_i (1/2 m_i v_i^2) return ((2 * self.ekin) / ( - 3 * dof * data.nist.BOLTZMANN / data.nist.HARTREE2J)) + dof * data.nist.BOLTZMANN / data.nist.HARTREE2J)) def __iter__(self): self._step = 0 @@ -359,8 +363,8 @@ def __next__(self): if self.incore_anyway: self.frames.append(current_frame) - if self.energy_output is not None: - self._write_energy() + if self.data_output is not None: + self._write_data() if self.trajectory_output is not None: self._write_coord() @@ -388,14 +392,15 @@ def _next(self): ''' raise NotImplementedError('Method Not Implemented') - def _write_energy(self): - '''Writes out the potential, kinetic, and total energy to the - self.energy_output stream. ''' + def _write_data(self): + '''Writes out the potential, kinetic, and total energy, temperature to the + self.data_output stream. ''' - output = '%8.2f %.12E %.12E %.12E' % (self.time, - self.epot, - self.ekin, - self.ekin + self.epot) + output = '%8.2f %.12E %.12E %.12E %3.4f' % (self.time, + self.epot, + self.ekin, + self.ekin + self.epot, + self.temperature()) # We follow OM of writing all the states at the end of the line if getattr(self.scanner.base, 'e_states', None) is not None: @@ -403,10 +408,10 @@ def _write_energy(self): for e in self.scanner.base.e_states: output += ' %.12E' % e - self.energy_output.write(output + '\n') + self.data_output.write(output + '\n') # If we don't flush, there is a possibility of losing data - self.energy_output.flush() + self.data_output.flush() def _write_coord(self): '''Writes out the current geometry to the self.trajectory_output @@ -494,3 +499,108 @@ def _next_velocity(self, next_accel): necessary equations of motion for the velocity is v(t_i+1) = v(t_i) + 0.5(a(t_i+1) + a(t_i))''' return self.veloc + 0.5 * self.dt * (self.accel + next_accel) + + +class NVTBerendson(_Integrator): + '''Berendsen (constant N, V, T) molecular dynamics + + Args: + method : lib.GradScanner or rhf.GradientsMixin instance, or + has nuc_grad_method method. + Method by which to compute the energy gradients and energies + in order to propagate the equations of motion. Realistically, + it can be any callable object such that it returns the energy + and potential energy gradient when given a mol. + + T : float + Target temperature for the NVT Ensemble. Given in K. + + taut : float + Time constant for Berendsen temperature coupling. + Given in atomic units. + + Attributes: + accel : ndarray + Current acceleration for the simulation. Values are given + in atomic units (Bohr/a.u.^2). Dimensions is (natm, 3) such as + + [[x1, y1, z1], + [x2, y2, z2], + [x3, y3, z3]] + ''' + + def __init__(self, method, T, taut, **kwargs): + self.T = T + self.taut = taut + self.accel = None + super().__init__(method, **kwargs) + + def _next(self): + '''Computes the next frame of the simulation and sets all internal + variables to this new frame. First computes the new geometry, + then the next acceleration, and finally the velocity, all according + to the Velocity Verlet algorithm. + + Returns: + The next frame of the simulation. + ''' + + # If no acceleration, compute that first, and then go + # onto the next step + if self.accel is None: + next_epot, next_accel = self._compute_accel() + + else: + self._scale_velocities() + self.mol.set_geom_(self._next_geometry(), unit='B') + self.mol.build() + next_epot, next_accel = self._compute_accel() + self.veloc = self._next_velocity(next_accel) + + self.epot = next_epot + self.ekin = self.compute_kinetic_energy() + self.accel = next_accel + + return _toframe(self) + + def _compute_accel(self): + '''Given the current geometry, computes the acceleration + for each atom.''' + e_tot, grad = self.scanner(self.mol) + if not self.scanner.converged: + raise RuntimeError('Gradients did not converge!') + + a = -1 * grad / self._masses.reshape(-1, 1) + return e_tot, a + + def _scale_velocities(self): + '''NVT Berendsen velocity scaling + v_rescale(t) = v(t) * (1 + ((T_target/T - 1) + * (/delta t / taut)))^(0.5) + ''' + tautscl = self.dt / self.taut + scl_temp = np.sqrt(1.0 + (self.T / self.temperature() - 1.0) * tautscl) + + # Limit the velocity scaling to reasonable values + # (taken from ase md/nvtberendson.py) + if scl_temp > 1.1: + scl_temp = 1.1 + if scl_temp < 0.9: + scl_temp = 0.9 + + self.veloc = self.veloc * scl_temp + return + + def _next_geometry(self): + '''Computes the next geometry using the Velocity Verlet algorithm. The + necessary equations of motion for the position is + r(t_i+1) = r(t_i) + /delta t * v(t_i) + 0.5(/delta t)^2 a(t_i) + ''' + return self.mol.atom_coords() + self.dt * self.veloc + \ + 0.5 * (self.dt ** 2) * self.accel + + def _next_velocity(self, next_accel): + '''Compute the next velocity using the Velocity Verlet algorithm. The + necessary equations of motion for the velocity is + v(t_i+1) = v(t_i) + /delta t * 0.5(a(t_i+1) + a(t_i))''' + return self.veloc + 0.5 * self.dt * (self.accel + next_accel) diff --git a/pyscf/md/test/test_NVTBerendson.py b/pyscf/md/test/test_NVTBerendson.py new file mode 100644 index 0000000000..4e6195d6c0 --- /dev/null +++ b/pyscf/md/test/test_NVTBerendson.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python +# Copyright 2014-2022 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Author: Aniruddha Seal + +import unittest +import numpy as np +from pyscf import gto, scf +import pyscf.md as md + +CHECK_STABILITY = False + +def setUpModule(): + global h2o, hf_scanner + h2o = gto.M(verbose=3, + output='/dev/null', + atom='''O -2.9103342153 -15.4805607073 -14.9344021104 + H -2.5833611256 -14.8540450112 -15.5615823519 + H -2.7404195919 -16.3470417109 -15.2830799053''', + basis='def2-svp') + + hf_scanner = scf.RHF(h2o) + hf_scanner.build() + hf_scanner.conv_tol_grad = 1e-6 + hf_scanner.max_cycle = 700 + +def tearDownModule(): + global h2o, hf_scanner + hf_scanner.stdout.close() + del h2o, hf_scanner + +class KnownValues(unittest.TestCase): + + def test_hf_water_init_veloc(self): + init_veloc = np.array([[-3.00670625e-05, -2.47219610e-04, -2.99235779e-04], + [-5.66022419e-05, -9.83256521e-04, -8.35299245e-04], + [-4.38260913e-04, -3.17970694e-04, 5.07818817e-04]]) + + driver = md.integrators.NVTBerendson(hf_scanner, veloc=init_veloc, + dt=20, steps=20, T=300, taut=413) + + driver.kernel() + self.assertAlmostEqual(driver.ekin, 4.244286252900E-03, 8) + self.assertAlmostEqual(driver.epot, -7.596117676337E+01, 8) + + final_coord = np.array([[ -5.51486188, -29.3425402, -28.32832762], + [ -4.8843336, -28.48585797, -29.75984939], + [ -5.30517791, -31.05471672, -28.76717135]]) + + self.assertTrue(np.allclose(driver.mol.atom_coords(), final_coord)) + + if CHECK_STABILITY: + + driver.steps = 990 + driver.kernel() + self.assertAlmostEqual(driver.temperature(), driver.T, delta=5) + +if __name__ == "__main__": + print("Full Tests for NVT Berendson Thermostat") + unittest.main() diff --git a/pyscf/md/test/test_prop.py b/pyscf/md/test/test_prop.py new file mode 100644 index 0000000000..964f7c698d --- /dev/null +++ b/pyscf/md/test/test_prop.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python +# Copyright 2014-2022 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Author: Aniruddha Seal + +import unittest +import numpy as np +from pyscf import gto, scf, data +import pyscf.md as md + +def setUpModule(): + global h2o, h2o_scanner, o2, o2_scanner + h2o = gto.M(verbose=3, + output='/dev/null', + atom='''O -2.9103342153 -15.4805607073 -14.9344021104 + H -2.5833611256 -14.8540450112 -15.5615823519 + H -2.7404195919 -16.3470417109 -15.2830799053''', + basis='def2-svp') + + h2o_scanner = scf.RHF(h2o) + h2o_scanner.build() + h2o_scanner.conv_tol_grad = 1e-6 + h2o_scanner.max_cycle = 700 + + o2 = gto.M(verbose=3, + output='/dev/null', + atom='''O 0 0 0.0298162549 + O 0 0 1.1701837504''', + basis='def2-svp') + + o2_scanner = scf.RHF(o2) + o2_scanner.build() + o2_scanner.conv_tol_grad = 1e-6 + o2_scanner.max_cycle = 700 + +def tearDownModule(): + global h2o, h2o_scanner, o2, o2_scanner + h2o_scanner.stdout.close() + o2_scanner.stdout.close() + del h2o, h2o_scanner, o2, o2_scanner + +class KnownValues(unittest.TestCase): + + def test_temperature_non_linear(self): + # Property Checked: Non-Linear Molecule Temperature + # unit-converted velocities temp = 298.13 obtained from ORCA + init_veloc = np.array([[-3.00670625e-05, -2.47219610e-04, -2.99235779e-04], + [-5.66022419e-05, -9.83256521e-04, -8.35299245e-04], + [-4.38260913e-04, -3.17970694e-04, 5.07818817e-04]]) + + driver = md.integrators.NVTBerendson(h2o_scanner, veloc=init_veloc, + dt=20, steps=50, T=300, taut=413) + + driver._masses = np.array( + [data.elements.COMMON_ISOTOPE_MASSES[m] * data.nist.AMU2AU + for m in driver.mol.atom_charges()]) + driver.ekin = driver.compute_kinetic_energy() + self.assertAlmostEqual(driver.temperature(), 298.13, delta=0.2) + + def test_temperature_linear(self): + # Property Checked: Linear Molecule Temperature + # unit-converted velocities temp = 468.31 obtained from ORCA + init_veloc = np.array([[-0., 0., 0.00039058], + [ 0., -0., -0.00039058]]) + + driver = md.integrators.NVTBerendson(o2_scanner, veloc=init_veloc, + dt=20, steps=50, T=300, taut=413) + + driver._masses = np.array( + [data.elements.COMMON_ISOTOPE_MASSES[m] * data.nist.AMU2AU + for m in driver.mol.atom_charges()]) + driver.ekin = driver.compute_kinetic_energy() + self.assertAlmostEqual(driver.temperature(), 468.31, delta=0.2) + +if __name__ == "__main__": + print("Property Tests") + unittest.main() From 96b2bf9365d0e2c762caaeac151433eba28fec15 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Fri, 15 Sep 2023 18:07:46 +0300 Subject: [PATCH 2/6] Fix lz_adaptation example; needs to run in D2h symmetry to work --- examples/symm/33-lz_adaption.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/symm/33-lz_adaption.py b/examples/symm/33-lz_adaption.py index 334a83b7c1..38589c5b3c 100644 --- a/examples/symm/33-lz_adaption.py +++ b/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() From 9e005f88d2d1ea79664266b63f5c36f432483e27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Chapoton?= Date: Fri, 15 Sep 2023 19:59:47 +0200 Subject: [PATCH 3/6] new collection of typos everywhere --- FEATURES | 2 +- examples/1-advanced/033-constrained_dft.py | 2 +- examples/dft/11-grid_scheme.py | 4 ++-- examples/local_orb/pmloc.py | 2 +- examples/scf/32-break_spin_symm.py | 2 +- pyscf/df/addons.py | 2 +- pyscf/df/incore.py | 2 +- pyscf/dft/libxc.py | 12 ++++++------ pyscf/dft/numint.py | 4 ++-- pyscf/dft/rks.py | 2 +- pyscf/eph/rhf.py | 2 +- pyscf/fci/cistring.py | 2 +- pyscf/fci/direct_spin1_cyl_sym.py | 2 +- pyscf/fci/direct_spin1_symm.py | 2 +- pyscf/fci/selected_ci.py | 2 +- pyscf/fci/selected_ci_symm.py | 2 +- pyscf/fci/spin_op.py | 2 +- pyscf/geomopt/berny_solver.py | 4 ++-- pyscf/geomopt/geometric_solver.py | 4 ++-- pyscf/grad/lagrange.py | 2 +- pyscf/gw/rpa.py | 4 ++-- pyscf/gw/urpa.py | 4 ++-- pyscf/lib/dft/libxc_itrf.c | 2 +- pyscf/lib/gto/nr_ecp.c | 2 +- pyscf/lib/pbc/ft_ao.c | 2 +- pyscf/mcscf/addons.py | 10 +++++----- pyscf/mcscf/casci.py | 2 +- pyscf/symm/addons.py | 4 ++-- pyscf/symm/basis.py | 2 +- pyscf/symm/geom.py | 11 ++++++----- 30 files changed, 51 insertions(+), 50 deletions(-) diff --git a/FEATURES b/FEATURES index 0d2caa2dbe..213967f26e 100644 --- a/FEATURES +++ b/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 diff --git a/examples/1-advanced/033-constrained_dft.py b/examples/1-advanced/033-constrained_dft.py index 067243f200..c40982c777 100644 --- a/examples/1-advanced/033-constrained_dft.py +++ b/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 diff --git a/examples/dft/11-grid_scheme.py b/examples/dft/11-grid_scheme.py index b58cf249a0..457fd38b75 100644 --- a/examples/dft/11-grid_scheme.py +++ b/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()) diff --git a/examples/local_orb/pmloc.py b/examples/local_orb/pmloc.py index 1e18154549..6a9626e796 100644 --- a/examples/local_orb/pmloc.py +++ b/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 diff --git a/examples/scf/32-break_spin_symm.py b/examples/scf/32-break_spin_symm.py index e934ee2943..2511f757f3 100644 --- a/examples/scf/32-break_spin_symm.py +++ b/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. diff --git a/pyscf/df/addons.py b/pyscf/df/addons.py index b042963b3a..033c68609b 100644 --- a/pyscf/df/addons.py +++ b/pyscf/df/addons.py @@ -30,7 +30,7 @@ ETB_BETA = getattr(__config__, 'df_addons_aug_dfbasis', 2.0) FIRST_ETB_ELEMENT = getattr(__config__, 'df_addons_aug_start_at', 36) # 'Rb' -# For code compatiblity in python-2 and python-3 +# For code compatibility in python-2 and python-3 if sys.version_info >= (3,): unicode = str diff --git a/pyscf/df/incore.py b/pyscf/df/incore.py index f89215b8c3..c9399a413e 100644 --- a/pyscf/df/incore.py +++ b/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 diff --git a/pyscf/dft/libxc.py b/pyscf/dft/libxc.py index de025b16bf..8543d12727 100644 --- a/pyscf/dft/libxc.py +++ b/pyscf/dft/libxc.py @@ -709,12 +709,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)) @@ -1061,7 +1061,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. @@ -1486,7 +1486,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))) diff --git a/pyscf/dft/numint.py b/pyscf/dft/numint.py index c4754a2a7d..0053b9979e 100644 --- a/pyscf/dft/numint.py +++ b/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. diff --git a/pyscf/dft/rks.py b/pyscf/dft/rks.py index 0dfceae566..2315c39a04 100644 --- a/pyscf/dft/rks.py +++ b/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] diff --git a/pyscf/eph/rhf.py b/pyscf/eph/rhf.py index da3eb81d67..c570fd23aa 100644 --- a/pyscf/eph/rhf.py +++ b/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) diff --git a/pyscf/fci/cistring.py b/pyscf/fci/cistring.py index 803cc4c1b6..014ed6e829 100644 --- a/pyscf/fci/cistring.py +++ b/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] diff --git a/pyscf/fci/direct_spin1_cyl_sym.py b/pyscf/fci/direct_spin1_cyl_sym.py index a666a393f2..a9f7a8bb2a 100644 --- a/pyscf/fci/direct_spin1_cyl_sym.py +++ b/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. ''' diff --git a/pyscf/fci/direct_spin1_symm.py b/pyscf/fci/direct_spin1_symm.py index c59b4f096f..41af31a34b 100644 --- a/pyscf/fci/direct_spin1_symm.py +++ b/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 diff --git a/pyscf/fci/selected_ci.py b/pyscf/fci/selected_ci.py index 14aae7861c..ae1dca17e1 100644 --- a/pyscf/fci/selected_ci.py +++ b/pyscf/fci/selected_ci.py @@ -900,7 +900,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''' diff --git a/pyscf/fci/selected_ci_symm.py b/pyscf/fci/selected_ci_symm.py index 6d3578f9d8..ad9565d017 100644 --- a/pyscf/fci/selected_ci_symm.py +++ b/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 diff --git a/pyscf/fci/spin_op.py b/pyscf/fci/spin_op.py index 973e74be04..6169f20f86 100644 --- a/pyscf/fci/spin_op.py +++ b/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 ''' diff --git a/pyscf/geomopt/berny_solver.py b/pyscf/geomopt/berny_solver.py index 4e039e12f8..0f477fc8e9 100644 --- a/pyscf/geomopt/berny_solver.py +++ b/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) diff --git a/pyscf/geomopt/geometric_solver.py b/pyscf/geomopt/geometric_solver.py index f11c46f9c6..0c39f9f0e9 100644 --- a/pyscf/geomopt/geometric_solver.py +++ b/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 diff --git a/pyscf/grad/lagrange.py b/pyscf/grad/lagrange.py index eb66360323..372ae9a949 100644 --- a/pyscf/grad/lagrange.py +++ b/pyscf/grad/lagrange.py @@ -114,7 +114,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, diff --git a/pyscf/gw/rpa.py b/pyscf/gw/rpa.py index 7a157f40bc..94f5e6d06e 100755 --- a/pyscf/gw/rpa.py +++ b/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) @@ -228,7 +228,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: diff --git a/pyscf/gw/urpa.py b/pyscf/gw/urpa.py index cc1324e777..5d338a7638 100755 --- a/pyscf/gw/urpa.py +++ b/pyscf/gw/urpa.py @@ -75,7 +75,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) @@ -159,7 +159,7 @@ def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40, x0=0.5): mo_energy : 2D array (2, nmo), mean-field mo energy mo_coeff : 3D array (2, nmo, nmo), mean-field mo coefficient Lpq : 4D array (2, naux, nmo, nmo), 3-index ERI - nw: interger, grid number + nw: integer, grid number x0: real, scaling factor for frequency grid Returns: diff --git a/pyscf/lib/dft/libxc_itrf.c b/pyscf/lib/dft/libxc_itrf.c index 989d4519d6..60641e5e43 100644 --- a/pyscf/lib/dft/libxc_itrf.c +++ b/pyscf/lib/dft/libxc_itrf.c @@ -356,7 +356,7 @@ static void _eval_xc(xc_func_type *func_x, int spin, int np, } break; default: - fprintf(stderr, "functional %d '%s' is not implmented\n", + fprintf(stderr, "functional %d '%s' is not implemented\n", func_x->info->number, func_x->info->name); raise_error; } diff --git a/pyscf/lib/gto/nr_ecp.c b/pyscf/lib/gto/nr_ecp.c index 5af611c83c..0385278f91 100644 --- a/pyscf/lib/gto/nr_ecp.c +++ b/pyscf/lib/gto/nr_ecp.c @@ -5320,7 +5320,7 @@ static int check_3c_overlap(int *shls, int *atm, int *bas, double *env, kprim = ecpbas[csh*BAS_SLOTS+NPRIM_OF]; ak = env + ecpbas[csh*BAS_SLOTS+PTR_EXP]; - // Test the last primitive funciton only because the basis + // Test the last primitive function only because the basis // functions are sorted so that the last one has the smallest // exponent. aijk = aij + ak[kprim-1]; diff --git a/pyscf/lib/pbc/ft_ao.c b/pyscf/lib/pbc/ft_ao.c index 4f38e522c1..720f0302bc 100644 --- a/pyscf/lib/pbc/ft_ao.c +++ b/pyscf/lib/pbc/ft_ao.c @@ -846,7 +846,7 @@ void PBCsupmol_ovlp_mask(int8_t *out, double cutoff, rr_cutoff = (log_a1 - log_cutoff) / a1; if (li > 0 && lj > 0 && a1 < 0.3) { // the contribution of r^n should be considered for - // overlap of smooth basis funcitons + // overlap of smooth basis functions li_a1 = li / -a1; lj_a1 = lj / -a1; for (i = i0; i < i1; i++) { diff --git a/pyscf/mcscf/addons.py b/pyscf/mcscf/addons.py index 1cd3cfa878..707236b623 100644 --- a/pyscf/mcscf/addons.py +++ b/pyscf/mcscf/addons.py @@ -377,7 +377,7 @@ def project_init_guess (casscf, mo_init, prev_mol=None, priority=None, use_hf_co Kwargs: prev_mol : an instance of :class:`Mole` - If given, the inital guess orbitals are associated to the + If given, the initial guess orbitals are associated to the basis of prev_mol. Otherwise, the orbitals are presumed to be in the basis of casscf.mol. Beware linear dependencies if you are projecting from a LARGER basis to a SMALLER one. @@ -555,7 +555,7 @@ def project_init_guess_old(casscf, init_mo, prev_mol=None): Kwargs: prev_mol : an instance of :class:`Mole` - If given, the inital guess orbitals are associated to the geometry + If given, the initial guess orbitals are associated to the geometry and basis of prev_mol. Otherwise, the orbitals are based of the geometry and basis of casscf.mol @@ -778,11 +778,11 @@ def get_fock(casscf, mo_coeff=None, ci=None): return casscf.get_fock(mo_coeff, ci) def cas_natorb(casscf, mo_coeff=None, ci=None, sort=False): - '''Natrual orbitals in CAS space + '''Natural orbitals in CAS space ''' if mo_coeff is None: mo_coeff = casscf.mo_coeff if _is_uhf_mo(mo_coeff): - raise RuntimeError('TODO: UCAS natrual orbitals') + raise RuntimeError('TODO: UCAS natural orbitals') else: return casscf.cas_natorb(mo_coeff, ci, sort=sort) @@ -851,7 +851,7 @@ class StateAverageMCSCFSolver: pass def state_average(casscf, weights=(0.5,0.5), wfnsym=None): - ''' State average over the energy. The energy funcitonal is + ''' State average over the energy. The energy functional is E = w1 + w2 + ... Note we may need change the FCI solver to diff --git a/pyscf/mcscf/casci.py b/pyscf/mcscf/casci.py index fee6c41ece..7b919ae70f 100644 --- a/pyscf/mcscf/casci.py +++ b/pyscf/mcscf/casci.py @@ -254,7 +254,7 @@ def get_fock(mc, mo_coeff=None, ci=None, eris=None, casdm1=None, verbose=None): def cas_natorb(mc, mo_coeff=None, ci=None, eris=None, sort=False, casdm1=None, verbose=None, with_meta_lowdin=WITH_META_LOWDIN): - '''Transform active orbitals to natrual orbitals, and update the CI wfn + '''Transform active orbitals to natural orbitals, and update the CI wfn accordingly Args: diff --git a/pyscf/symm/addons.py b/pyscf/symm/addons.py index 22fdbf8c27..edb58f1c21 100644 --- a/pyscf/symm/addons.py +++ b/pyscf/symm/addons.py @@ -310,13 +310,13 @@ def symmetrize_multidim(mol, mo, s=None, keep_phase=getattr(__config__, 'symm_addons_symmetrize_multidim_keep_phase', True)): '''Symmetrize orbitals with respect to multidimensional irreps. - Make coefficients of partner funtions of multidimensional irreps to be the same. + Make coefficients of partner functions of multidimensional irreps to be the same. The functions uses the convention of the libmsym interface, that introduces underscores to the labels of multidimensional irreps partners. Args: mol : an instance of :class:`Mole` - Symmerty-adapted basis with multidimensional irreps should be generated by libmsym + Symmetry-adapted basis with multidimensional irreps should be generated by libmsym mo : 2D float array The orbital space to symmetrize diff --git a/pyscf/symm/basis.py b/pyscf/symm/basis.py index 982bc54367..321fa5bbf9 100644 --- a/pyscf/symm/basis.py +++ b/pyscf/symm/basis.py @@ -209,7 +209,7 @@ def _basis_offset_for_atoms(atoms, basis_tab): def _num_contract(basis): if isinstance(basis[1], int): - # This branch should never be reached if basis_tab is formated by + # This branch should never be reached if basis_tab is formatted by # function mole.format_basis nctr = len(basis[2]) - 1 else: diff --git a/pyscf/symm/geom.py b/pyscf/symm/geom.py index 2ee075e63a..1a81ff820d 100644 --- a/pyscf/symm/geom.py +++ b/pyscf/symm/geom.py @@ -49,10 +49,11 @@ TOLERANCE = getattr(__config__, 'symm_geom_tol', 1e-5) -# For code compatiblity in python-2 and python-3 +# For code compatibility in python-2 and python-3 if sys.version_info >= (3,): unicode = str + def parallel_vectors(v1, v2, tol=TOLERANCE): if numpy.allclose(v1, 0, atol=tol) or numpy.allclose(v2, 0, atol=tol): return True @@ -121,7 +122,7 @@ def alias_axes(axes, ref): def _adjust_planar_c2v(atom_coords, axes): '''Adjust axes for planar molecules''' # Following http://iopenshell.usc.edu/resources/howto/symmetry/ - # See also dicussions in issue #1201 + # See also discussions in issue #1201 # * planar C2v molecules should be oriented such that the X axis is perpendicular # to the plane of the molecule, and the Z axis is the axis of symmetry; natm = len(atom_coords) @@ -135,7 +136,7 @@ def _adjust_planar_c2v(atom_coords, axes): def _adjust_planar_d2h(atom_coords, axes): '''Adjust axes for planar molecules''' # Following http://iopenshell.usc.edu/resources/howto/symmetry/ - # See also dicussions in issue #1201 + # See also discussions in issue #1201 # * planar D2h molecules should be oriented such that the X axis is # perpendicular to the plane of the molecule, and the Z axis passes through # the greatest number of atoms. @@ -553,7 +554,7 @@ class RotationAxisNotFound(PointGroupSymmetryError): class SymmSys(object): def __init__(self, atoms, basis=None): self.atomtypes = mole.atom_types(atoms, basis) - # fake systems, which treates the atoms of different basis as different atoms. + # fake systems, which treats the atoms of different basis as different atoms. # the fake systems do not have the same symmetry as the potential # it's only used to determine the main (Z-)axis chg1 = numpy.pi - 2 @@ -614,7 +615,7 @@ def symmetric_for(self, op): for lst in self.group_atoms_by_distance: r0 = self.atoms[lst,1:] r1 = numpy.dot(r0, op) -# FIXME: compare whehter two sets of coordinates are identical + # FIXME: compare whether two sets of coordinates are identical yield all((_vec_in_vecs(x, r0) for x in r1)) def has_icenter(self): From 0449778dfb515fb0a240e2fd8f308a0b3eaa697c Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Fri, 15 Sep 2023 17:05:33 -0700 Subject: [PATCH 4/6] Optimize grad_nuc --- pyscf/grad/rhf.py | 20 +++++++++----------- pyscf/grad/test/test_rhf.py | 19 ++++++++++++++++++- 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/pyscf/grad/rhf.py b/pyscf/grad/rhf.py index 87ef29c468..6d95dc6de2 100644 --- a/pyscf/grad/rhf.py +++ b/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) diff --git a/pyscf/grad/test/test_rhf.py b/pyscf/grad/test/test_rhf.py index 78342aadb4..108d4a56e2 100644 --- a/pyscf/grad/test/test_rhf.py +++ b/pyscf/grad/test/test_rhf.py @@ -201,8 +201,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() - From 868dac2857a022a3f718cf0fc6598ccc477440e8 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Sat, 16 Sep 2023 16:26:55 -0500 Subject: [PATCH 5/6] Newton casscf hcc (issue #701) (#1869) * Hessian_CC in newton_casscf (issue #731) Hessian-vector product should have zero overlap w/ keyframe CI vector at all times. Remove factors of 2 which spoiled this. * update newton_casscf unit tests --------- Co-authored-by: Matthew R Hermes --- pyscf/fci/rdm.py | 4 ++-- pyscf/mcscf/newton_casscf.py | 8 +++++--- pyscf/mcscf/test/test_newton_casscf.py | 8 ++++---- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/pyscf/fci/rdm.py b/pyscf/fci/rdm.py index 58be01c4fe..284053c014 100644 --- a/pyscf/fci/rdm.py +++ b/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), diff --git a/pyscf/mcscf/newton_casscf.py b/pyscf/mcscf/newton_casscf.py index 941def1e56..b6e4d9c352 100644 --- a/pyscf/mcscf/newton_casscf.py +++ b/pyscf/mcscf/newton_casscf.py @@ -293,7 +293,8 @@ def g_update(u, fcivec): # h_diag[ncore:nocc,ncore:nocc] -= v_diag[:,ncore:nocc]*2 h_diag = casscf.pack_uniq_var(h_diag) - hci_diag = [hd - ec - gc*c*4 for hd, ec, gc, c in zip (_Hdiag (h1cas_0, eri_cas), eci0, gci, ci0)] + # Fix intermediate normalization? - MRH 2023/09/15 + hci_diag = [hd - ec - gc*c*2 for hd, ec, gc, c in zip (_Hdiag (h1cas_0, eri_cas), eci0, gci, ci0)] hci_diag = [h * w for h, w in zip (hci_diag, weights)] hdiag_all = numpy.hstack((h_diag*2, _pack_ci (hci_diag)*2)) @@ -308,9 +309,10 @@ def h_op(x): ci1 = _unpack_ci (x[ngorb:]) # H_cc + # Fix intermediate normalization? - MRH 2023/09/15 hci1 = [hc1 - c1 * ec0 for hc1, c1, ec0 in zip (_Hci (h1cas_0, eri_cas, ci1), ci1, eci0)] - hci1 = [hc1 - 2*(hc0 - c0*ec0)*c0.dot(c1) for hc1, hc0, c0, ec0, c1 in zip (hci1, hci0, ci0, eci0, ci1)] - hci1 = [hc1 - 2*c0*(hc0 - c0*ec0).dot(c1) for hc1, hc0, c0, ec0, c1 in zip (hci1, hci0, ci0, eci0, ci1)] + hci1 = [hc1 - (hc0 - c0*ec0)*c0.dot(c1) for hc1, hc0, c0, ec0, c1 in zip (hci1, hci0, ci0, eci0, ci1)] + hci1 = [hc1 - c0*(hc0 - c0*ec0).dot(c1) for hc1, hc0, c0, ec0, c1 in zip (hci1, hci0, ci0, eci0, ci1)] # H_co rc = x1[:,:ncore] diff --git a/pyscf/mcscf/test/test_newton_casscf.py b/pyscf/mcscf/test/test_newton_casscf.py index db4e7c15a7..f2d94c3e9a 100644 --- a/pyscf/mcscf/test/test_newton_casscf.py +++ b/pyscf/mcscf/test/test_newton_casscf.py @@ -88,11 +88,11 @@ def test_gen_g_hop(self): ci0/= numpy.linalg.norm(ci0) gall, gop, hop, hdiag = newton_casscf.gen_g_hop(mc, mo, ci0, mc.ao2mo(mo)) self.assertAlmostEqual(lib.fp(gall), 21.288022525148595, 8) - self.assertAlmostEqual(lib.fp(hdiag), -4.6864640132374618, 8) + self.assertAlmostEqual(lib.fp(hdiag), -15.618395788969842, 8) x = numpy.random.random(gall.size) u, ci1 = newton_casscf.extract_rotation(mc, x, 1, ci0) self.assertAlmostEqual(lib.fp(gop(u, ci1)), -412.9441873541524, 8) - self.assertAlmostEqual(lib.fp(hop(x)), 73.358310983341198, 8) + self.assertAlmostEqual(lib.fp(hop(x)), 24.04569925660985, 8) def test_get_grad(self): # mc.e_tot may be converged to -3.6268060853430573 @@ -107,11 +107,11 @@ def test_sa_gen_g_hop(self): ci0 = list (ci0.reshape ((2,6,6))) gall, gop, hop, hdiag = newton_casscf.gen_g_hop(sa, mo, ci0, sa.ao2mo(mo)) self.assertAlmostEqual(lib.fp(gall), 32.46973284682045, 8) - self.assertAlmostEqual(lib.fp(hdiag), -63.6527761153809, 8) + self.assertAlmostEqual(lib.fp(hdiag), -70.61862254321517, 8) x = numpy.random.random(gall.size) u, ci1 = newton_casscf.extract_rotation(sa, x, 1, ci0) self.assertAlmostEqual(lib.fp(gop(u, ci1)), -49.017079186126, 8) - self.assertAlmostEqual(lib.fp(hop(x)), 169.47893548740288, 8) + self.assertAlmostEqual(lib.fp(hop(x)), 136.2077988624156, 8) def test_sa_get_grad(self): self.assertAlmostEqual(sa.e_tot, -3.62638372957158, 7) From f3334e68fc1753aac548ad0c2af42eed370346a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Chapoton?= Date: Mon, 14 Aug 2023 14:19:48 +0200 Subject: [PATCH 6/6] fix other rst issues --- pyscf/pbc/ao2mo/eris.py | 15 +++++++++++---- pyscf/pbc/gto/cell.py | 1 + pyscf/pbc/gto/eval_gto.py | 2 +- pyscf/pbc/scf/hf.py | 4 ++-- pyscf/pbc/symm/space_group.py | 2 ++ pyscf/pbc/symm/symmetry.py | 4 ++++ pyscf/pbc/tools/pbc.py | 12 +++++++----- pyscf/pbc/tools/pywannier90.py | 28 ++++++++++++++++++++++------ pyscf/scf/addons.py | 3 ++- pyscf/sgx/sgx_jk.py | 6 ++++-- pyscf/symm/sph.py | 4 +++- pyscf/tools/fcidump.py | 7 +++++-- 12 files changed, 64 insertions(+), 24 deletions(-) diff --git a/pyscf/pbc/ao2mo/eris.py b/pyscf/pbc/ao2mo/eris.py index 0eb011b085..a21077bb1c 100644 --- a/pyscf/pbc/ao2mo/eris.py +++ b/pyscf/pbc/ao2mo/eris.py @@ -58,13 +58,16 @@ def get_mo_eri(cell, mo_coeffs, kpts=None): def get_mo_pairs_G(cell, mo_coeffs, kpts=None, q=None): '''Calculate forward (G|ij) FFT of all MO pairs. + TODO: - Implement simplifications for real orbitals. + Args: mo_coeff: length-2 list of (nao,nmo) ndarrays The two sets of MO coefficients to use in calculating the - product |ij). + product ``|ij)``. + Returns: - mo_pairs_G : (ngrids, nmoi*nmoj) ndarray + mo_pairs_G : ``(ngrids, nmoi*nmoj)`` ndarray The FFT of the real-space MO pairs. ''' coords = gen_uniform_grids(cell) @@ -108,13 +111,16 @@ def get_mo_pairs_G(cell, mo_coeffs, kpts=None, q=None): def get_mo_pairs_invG(cell, mo_coeffs, kpts=None, q=None): '''Calculate "inverse" (ij|G) FFT of all MO pairs. + TODO: - Implement simplifications for real orbitals. + Args: mo_coeff: length-2 list of (nao,nmo) ndarrays The two sets of MO coefficients to use in calculating the - product |ij). + product ``|ij)``. + Returns: - mo_pairs_invG : (ngrids, nmoi*nmoj) ndarray + mo_pairs_invG : ``(ngrids, nmoi*nmoj)`` ndarray The inverse FFTs of the real-space MO pairs. ''' coords = gen_uniform_grids(cell) @@ -158,6 +164,7 @@ def get_mo_pairs_invG(cell, mo_coeffs, kpts=None, q=None): def assemble_eri(cell, orb_pair_invG1, orb_pair_G2, q=None, verbose=logger.INFO): '''Assemble 4-index electron repulsion integrals. + Returns: (nmo1*nmo2, nmo3*nmo4) ndarray ''' diff --git a/pyscf/pbc/gto/cell.py b/pyscf/pbc/gto/cell.py index 8fc38d5e61..1cd47279fc 100644 --- a/pyscf/pbc/gto/cell.py +++ b/pyscf/pbc/gto/cell.py @@ -1689,6 +1689,7 @@ def get_scaled_kpts(self, abs_kpts, kpts_in_ibz=True): Args: abs_kpts : (nkpts, 3) ndarray of floats or :class:`KPoints` object + kpts_in_ibz : bool If True, return k-points in IBZ; otherwise, return k-points in BZ. Default value is True. This has effects only if abs_kpts is a diff --git a/pyscf/pbc/gto/eval_gto.py b/pyscf/pbc/gto/eval_gto.py index 0c598eba5a..723c389209 100644 --- a/pyscf/pbc/gto/eval_gto.py +++ b/pyscf/pbc/gto/eval_gto.py @@ -34,7 +34,7 @@ def eval_gto(cell, eval_name, coords, comp=None, kpts=None, kpt=None, r'''Evaluate PBC-AO function value on the given grids, Args: - eval_name : str + eval_name : str:: ========================== ======================= Function Expression diff --git a/pyscf/pbc/scf/hf.py b/pyscf/pbc/scf/hf.py index 3a9c2b2161..85c98dfe97 100644 --- a/pyscf/pbc/scf/hf.py +++ b/pyscf/pbc/scf/hf.py @@ -120,7 +120,7 @@ def get_j(cell, dm, hermi=1, vhfopt=None, kpt=np.zeros(3), kpts_band=None): kpt : (3,) ndarray The "inner" dummy k-point at which the DM was evaluated (or sampled). - kpts_band : (3,) ndarray or (*,3) ndarray + kpts_band : ``(3,)`` ndarray or ``(*,3)`` ndarray An arbitrary "band" k-point at which J is evaluated. Returns: @@ -150,7 +150,7 @@ def get_jk(mf, cell, dm, hermi=1, vhfopt=None, kpt=np.zeros(3), kpt : (3,) ndarray The "inner" dummy k-point at which the DM was evaluated (or sampled). - kpts_band : (3,) ndarray or (*,3) ndarray + kpts_band : ``(3,)`` ndarray or ``(*,3)`` ndarray An arbitrary "band" k-point at which J and K are evaluated. Returns: diff --git a/pyscf/pbc/symm/space_group.py b/pyscf/pbc/symm/space_group.py index e70030a402..8fdfcbc588 100644 --- a/pyscf/pbc/symm/space_group.py +++ b/pyscf/pbc/symm/space_group.py @@ -250,8 +250,10 @@ def __str__(self): class SpaceGroup(lib.StreamObject): ''' Determines the space group of a lattice. + Attributes: cell : :class:`Cell` object + symprec : float Numerical tolerance for determining the space group. Default value is 1e-6 in the unit of length. diff --git a/pyscf/pbc/symm/symmetry.py b/pyscf/pbc/symm/symmetry.py index 3618039c8f..c79bc81167 100644 --- a/pyscf/pbc/symm/symmetry.py +++ b/pyscf/pbc/symm/symmetry.py @@ -32,6 +32,7 @@ def get_Dmat(op, l): ''' Get Wigner D-matrix + Args: op : (3,3) ndarray rotation operator in (x,y,z) system @@ -133,7 +134,9 @@ class Symmetry(): Attributes: cell : :class:`Cell` object + spacegroup : :class:`SpaceGroup` object + symmorphic : bool Whether space group is symmorphic has_inversion : bool @@ -288,6 +291,7 @@ def transform_mo_coeff(cell, kpt_scaled, mo_coeff, op, Dmats): Args: cell : :class:`Cell` object + kpt_scaled : (3,) array scaled k-point mo_coeff : (nao, nmo) array diff --git a/pyscf/pbc/tools/pbc.py b/pyscf/pbc/tools/pbc.py index 447f752590..3625ec5ff3 100644 --- a/pyscf/pbc/tools/pbc.py +++ b/pyscf/pbc/tools/pbc.py @@ -197,7 +197,7 @@ def get_coulG(cell, k=np.zeros(3), exx=False, mf=None, mesh=None, Gv=None, k : (3,) ndarray k-point exx : bool or str - Whether this is an exchange matrix element. + Whether this is an exchange matrix element mf : instance of :class:`SCF` Returns: @@ -206,8 +206,8 @@ def get_coulG(cell, k=np.zeros(3), exx=False, mf=None, mesh=None, Gv=None, mesh : (3,) ndarray of ints (= nx,ny,nz) The number G-vectors along each direction. omega : float - Enable Coulomb kernel erf(|omega|*r12)/r12 if omega > 0 - and erfc(|omega|*r12)/r12 if omega < 0. + Enable Coulomb kernel ``erf(|omega|*r12)/r12`` if omega > 0 + and ``erfc(|omega|*r12)/r12`` if omega < 0. Note this parameter is slightly different to setting cell.omega for the treatment of exxdiv (at G0). cell.omega affects Ewald probe charge at G0. It is used mostly by RSH functionals for @@ -553,12 +553,14 @@ def find_boundary(a): def super_cell(cell, ncopy, wrap_around=False): '''Create an ncopy[0] x ncopy[1] x ncopy[2] supercell of the input cell - Note this function differs from :fun:`cell_plus_imgs` that cell_plus_imgs + Note this function differs from :func:`cell_plus_imgs` that cell_plus_imgs creates images in both +/- direction. Args: cell : instance of :class:`Cell` + ncopy : (3,) array + wrap_around : bool Put the original cell centered on the super cell. It has the effects corresponding to the parameter wrap_around of @@ -596,7 +598,7 @@ def super_cell(cell, ncopy, wrap_around=False): def cell_plus_imgs(cell, nimgs): '''Create a supercell via nimgs[i] in each +/- direction, as in get_lattice_Ls(). - Note this function differs from :fun:`super_cell` that super_cell only + Note this function differs from :func:`super_cell` that super_cell only stacks the images in + direction. Args: diff --git a/pyscf/pbc/tools/pywannier90.py b/pyscf/pbc/tools/pywannier90.py index 0b1877af05..775a70faaa 100644 --- a/pyscf/pbc/tools/pywannier90.py +++ b/pyscf/pbc/tools/pywannier90.py @@ -22,8 +22,10 @@ email: pqh3.14@gmail.com (2) -Another wannier90 python interface is available on the repo: +Another wannier90 python interface is available on the repo:: + https://github.com/zhcui/wannier90 + Contact its author "Zhihao Cui" for more details of installation and implementations. ''' @@ -152,6 +154,7 @@ def R_r(r_norm, r=1, zona=1): def theta(func, cost, phi): r''' Basic angular functions (s,p,d,f) used to compute \Theta_{l,m_r}(\theta,\phi) + ref: Table 3.1 of the Wannier90 User guide Link: https://github.com/wannier-developers/wannier90/raw/v3.1.0/doc/compiled_docs/user_guide.pdf ''' @@ -194,6 +197,7 @@ def theta(func, cost, phi): def theta_lmr(l, mr, cost, phi): r''' Compute the value of \Theta_{l,m_r}(\theta,\phi) + ref: Table 3.1 and 3.2 of the Wannier90 User guide Link: https://github.com/wannier-developers/wannier90/raw/v3.1.0/doc/compiled_docs/user_guide.pdf ''' @@ -282,11 +286,14 @@ def theta_lmr(l, mr, cost, phi): def g_r(grids_coor, site, l, mr, r, zona, x_axis=[1,0,0], z_axis=[0,0,1], unit='B'): r''' Evaluate the projection function g(r) or \Theta_{l,m_r}(\theta,\phi) on a grid + ref: Chapter 3, wannier90 User Guide + Attributes: grids_coor : a grids for the cell of interest site : absolute coordinate (in Borh/Angstrom) of the g(r) in the cell l, mr : l and mr value in the Table 3.1 and 3.2 of the ref + Return: theta_lmr : an array (ngrid, value) of g(r) @@ -317,6 +324,7 @@ def g_r(grids_coor, site, l, mr, r, zona, x_axis=[1,0,0], z_axis=[0,0,1], unit=' def get_wigner_seitz_supercell(w90, ws_search_size=[2,2,2], ws_distance_tol=1e-6): ''' Return a grid that contains all the lattice within the Wigner-Seitz supercell + Ref: the hamiltonian_wigner_seitz(count_pts) in wannier90/src/hamittonian.F90 ''' @@ -365,6 +373,7 @@ def get_wigner_seitz_supercell(w90, ws_search_size=[2,2,2], ws_distance_tol=1e-6 def R_wz_sc(w90, R_in, R0, ws_search_size=[2,2,2], ws_distance_tol=1e-6): ''' TODO: document it + Ref: This is the replication of the R_wz_sc function of ws_distance.F90 ''' ndegenx = 8 #max number of unit cells that can touch in a single point (i.e. vertex of cube) @@ -421,6 +430,7 @@ def R_wz_sc(w90, R_in, R0, ws_search_size=[2,2,2], ws_distance_tol=1e-6): def ws_translate_dist(w90, irvec, ws_search_size=[2,2,2], ws_distance_tol=1e-6): ''' TODO: document it + Ref: This is the replication of the ws_translate_dist function of ws_distance.F90 ''' nrpts = irvec.shape[0] @@ -565,7 +575,7 @@ def kernel(self, external_AME=None): def make_win(self): ''' - Make a basic *.win file for wannier90 + Make a basic ``*.win`` file for wannier90 ''' win_file = open('wannier90.win', "w") @@ -846,8 +856,9 @@ def interpolate_ham_kpts(self, frac_kpts, ham_kpts=None, use_ws_distance=True, ws_search_size=[2,2,2], ws_distance_tol=1e-6): ''' Interpolate the band structure using the Slater-Koster scheme - Return: - eigenvalues and eigenvectors at the desired kpts + + Return: + eigenvalues and eigenvectors at the desired kpts ''' assert self.U_matrix is not None, "You must wannierize first, then you can run this function" @@ -871,8 +882,9 @@ def interpolate_ham_kpts(self, frac_kpts, ham_kpts=None, def interpolate_band(self, frac_kpts, ham_kpts=None, use_ws_distance=True, ws_search_size=[2,2,2], ws_distance_tol=1e-6): ''' Interpolate the band structure using the Slater-Koster scheme - Return: - eigenvalues and eigenvectors at the desired kpts + + Return: + eigenvalues and eigenvectors at the desired kpts ''' assert self.U_matrix is not None, ( @@ -1031,7 +1043,9 @@ def get_guess_orb(self, frac_site=[0,0,0], l=0, mr=1, r=1, def plot_wf(self, outfile='MLWF', wf_list=None, supercell=[1,1,1], grid=[50,50,50]): ''' Export Wannier function at cell R + xsf format: http://web.mit.edu/xcrysden_v1.5.60/www/XCRYSDEN/doc/XSF.html + Attributes: wf_list : a list of MLWFs to plot supercell : a supercell used for plotting @@ -1090,7 +1104,9 @@ def plot_guess_orbs(self, outfile='guess_orb', frac_site=[0,0,0], l=0, mr=1, r=1 supercell=[1,1,1], grid=[50,50,50]): ''' Export Wannier function at cell R + xsf format: http://web.mit.edu/xcrysden_v1.5.60/www/XCRYSDEN/doc/XSF.html + Attributes: wf_list : a list of MLWFs to plot supercell : a supercell used for plotting diff --git a/pyscf/scf/addons.py b/pyscf/scf/addons.py index 986bf9f6b4..df0dc91311 100644 --- a/pyscf/scf/addons.py +++ b/pyscf/scf/addons.py @@ -918,8 +918,9 @@ def fast_newton(mf, mo_coeff=None, mo_occ=None, dm0=None, function first setup the initial guess from density fitting calculation then use for Newton solver and call Newton solver. + Newton solver attributes [max_cycle_inner, max_stepsize, ah_start_tol, - ah_conv_tol, ah_grad_trust_region, ...] can be passed through **newton_kwargs. + ah_conv_tol, ah_grad_trust_region, ...] can be passed through ``**newton_kwargs``. ''' import copy from pyscf.lib import logger diff --git a/pyscf/sgx/sgx_jk.py b/pyscf/sgx/sgx_jk.py index 31a5a83a62..b92f9995c9 100644 --- a/pyscf/sgx/sgx_jk.py +++ b/pyscf/sgx/sgx_jk.py @@ -252,8 +252,10 @@ def batch_nuc(mol, grid_coords, out=None): def _gen_jk_direct(mol, aosym, with_j, with_k, direct_scf_tol, sgxopt=None, pjs=False): '''Contraction between sgX Coulomb integrals and density matrices - J: einsum('guv,xg->xuv', gbn, dms) if dms == rho at grid - einsum('gij,xij->xg', gbn, dms) if dms are density matrices + + J: einsum('guv,xg->xuv', gbn, dms) if dms == rho at grid, + or einsum('gij,xij->xg', gbn, dms) if dms are density matrices + K: einsum('gtv,xgt->xgv', gbn, fg) ''' if sgxopt is None: diff --git a/pyscf/symm/sph.py b/pyscf/symm/sph.py index fe27623e01..59245de011 100644 --- a/pyscf/symm/sph.py +++ b/pyscf/symm/sph.py @@ -106,8 +106,10 @@ def multipoles(r, lmax, reorder_dipole=True): def sph_pure2real(l, reorder_p=True): r''' Transformation matrix: from the pure spherical harmonic functions Y_m to - the real spherical harmonic functions O_m. + the real spherical harmonic functions O_m:: + O_m = \sum Y_m' * U(m',m) + Y(-1) = 1/\sqrt(2){-iO(-1) + O(1)}; Y(1) = 1/\sqrt(2){-iO(-1) - O(1)} Y(-2) = 1/\sqrt(2){-iO(-2) + O(2)}; Y(2) = 1/\sqrt(2){iO(-2) + O(2)} O(-1) = i/\sqrt(2){Y(-1) + Y(1)}; O(1) = 1/\sqrt(2){Y(-1) - Y(1)} diff --git a/pyscf/tools/fcidump.py b/pyscf/tools/fcidump.py index 96b43fddf6..ca50a6a3ce 100644 --- a/pyscf/tools/fcidump.py +++ b/pyscf/tools/fcidump.py @@ -218,10 +218,13 @@ def read(filename, molpro_orbsym=MOLPRO_ORBSYM, verbose=True): Kwargs: molpro_orbsym (bool): Whether the orbsym in the FCIDUMP file is in - Molpro orbsym convention as documented in - https://www.molpro.net/info/current/doc/manual/node36.html + Molpro orbsym convention as documented in:: + + https://www.molpro.net/info/current/doc/manual/node36.html + In return, orbsym is converted to pyscf symmetry convention verbose (bool): Whether to print debugging information + ''' if verbose: print('Parsing %s' % filename)