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/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() 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()