Merge branch 'master' into newton_casscf_Hcc

CHANGELOG

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

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/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/md/__init__.py

@@ -45,4 +45,3 @@ def set_seed(seed):
 from pyscf.md import integrators, distributions
 
 NVE = integrators.VelocityVerlet
-

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 <[email protected]>
+# Authors: Matthew Hennefarth <[email protected]>,
+#          Aniruddha Seal <[email protected]>
 
 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,25 +392,26 @@ 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:
             if len(self.scanner.base.e_states) > 1:
                 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)

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 <[email protected]>
+
+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()