Commit d339f840 authored by Lauri Himanen's avatar Lauri Himanen

Added MD test, now parsing MD energy and temperature + statistics.

parent 4ab749dd
......@@ -377,15 +377,18 @@ class TestDFTGaussianMD(unittest.TestCase):
result = self.results["sampling_method"]
self.assertEqual(result, "molecular_dynamics")
# def test_number_of_atoms(self):
# result = self.results["number_of_atoms"]
# expected_result = np.array(5*[2])
# self.assertTrue(np.array_equal(result, expected_result))
def test_ensemble_type(self):
result = self.results["ensemble_type"]
self.assertEqual(result, "NVT")
def test_frame_sequence_to_sampling_ref(self):
result = self.results["frame_sequence_to_sampling_ref"]
self.assertEqual(result, 0)
def test_frame_sequence_local_frames_ref(self):
result = self.results["frame_sequence_local_frames_ref"]
self.assertTrue(np.array_equal(result, np.array(range(5))))
def test_number_of_frames_in_sequence(self):
result = self.results["number_of_frames_in_sequence"]
self.assertEqual(result, 5)
......@@ -398,31 +401,31 @@ class TestDFTGaussianMD(unittest.TestCase):
result = self.results["frame_sequence_kinetic_energy"]
self.assertTrue(np.array_equal(result, self.kin))
def test_frame_sequence_local_frames_ref(self):
result = self.results["frame_sequence_local_frames_ref"]
self.assertTrue(np.array_equal(result, np.array(range(5))))
def test_frame_sequence_temperature(self):
result = self.results["frame_sequence_temperature"]
self.assertTrue(np.array_equal(result, self.temp))
# def test_frame_sequence_conserved_quantity(self):
# result = self.results["frame_sequence_conserved_quantity"]
# self.assertTrue(np.array_equal(result, self.cons))
def test_frame_sequence_time(self):
result = self.results["frame_sequence_time"]
self.assertTrue(np.array_equal(result, self.time))
# def test_frame_sequence_temperature(self):
# result = self.results["frame_sequence_temperature"]
# self.assertTrue(np.array_equal(result, self.temp))
def test_frame_sequence_potential_energy_stats(self):
result = self.results["frame_sequence_potential_energy_stats"]
expected_result = np.array([self.pot.mean(), self.pot.std()])
self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.00001e-18))
self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.00001e-20))
# def test_frame_sequence_time(self):
# result = self.results["frame_sequence_time"]
# expected_result = convert_unit(
# np.array([
# 4,
# 8,
# 12,
# 16,
# 20,
# ]),
# "hbar/hartree"
# )
# self.assertTrue(np.array_equal(result, expected_result))
def test_frame_sequence_kinetic_energy_stats(self):
result = self.results["frame_sequence_kinetic_energy_stats"]
expected_result = np.array([self.kin.mean(), self.kin.std()])
self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.00001e-18))
self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.00001e-20))
def test_frame_sequence_temperature_stats(self):
result = self.results["frame_sequence_temperature_stats"]
expected_result = np.array([self.temp.mean(), self.temp.std()])
self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.001))
self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.0001))
# def test_atom_positions(self):
# result = self.results["atom_positions"]
......@@ -483,30 +486,6 @@ class TestDFTGaussianMD(unittest.TestCase):
# self.assertTrue(np.array_equal(result[0, :], expected_start))
# self.assertTrue(np.array_equal(result[-1, :], expected_end))
# def test_frame_sequence_potential_energy_stats(self):
# result = self.results["frame_sequence_potential_energy_stats"]
# expected_result = np.array([self.pot.mean(), self.pot.std()])
# self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.00001e-18))
# self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.00001e-20))
# def test_frame_sequence_kinetic_energy_stats(self):
# result = self.results["frame_sequence_kinetic_energy_stats"]
# expected_result = np.array([self.kin.mean(), self.kin.std()])
# self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.00001e-18))
# self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.00001e-20))
# def test_frame_sequence_conserved_quantity_stats(self):
# result = self.results["frame_sequence_conserved_quantity_stats"]
# expected_result = np.array([self.cons.mean(), self.cons.std()])
# self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.00001e-18))
# self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.00001e-20))
# def test_frame_sequence_temperature_stats(self):
# result = self.results["frame_sequence_temperature_stats"]
# expected_result = np.array([self.temp.mean(), self.temp.std()])
# self.assertTrue(np.allclose(result[0], expected_result[0], rtol=0, atol=0.001))
# self.assertTrue(np.allclose(result[1], expected_result[1], rtol=0, atol=0.0001))
#===============================================================================
class TestDFTGaussianXCFunctional(unittest.TestCase):
......@@ -571,7 +550,7 @@ if __name__ == '__main__':
suites = []
# suites.append(unittest.TestLoader().loadTestsFromTestCase(TestDFTGaussianEnergy))
# suites.append(unittest.TestLoader().loadTestsFromTestCase(TestDFTGaussianForce))
# suites.append(unittest.TestLoader().loadTestsFromTestCase(TestDFTGaussianGeoOpt))
suites.append(unittest.TestLoader().loadTestsFromTestCase(TestDFTGaussianGeoOpt))
# suites.append(unittest.TestLoader().loadTestsFromTestCase(TestDFTGaussianXCFunctional))
suites.append(unittest.TestLoader().loadTestsFromTestCase(TestDFTGaussianMD))
......
......@@ -19,7 +19,6 @@ class NWChemMainParser(MainHierarchicalParser):
super(NWChemMainParser, self).__init__(file_path, parser_context)
self.n_scf_iterations = 0
self.latest_dft_section = None
# self.frame_sequence_local_frames_ref = []
self.method_index = None
self.system_index = None
self.save_method = False
......@@ -39,6 +38,9 @@ class NWChemMainParser(MainHierarchicalParser):
self.frame_sequence_cache.add("frame_sequence_local_frames_ref", [], single=False, update=True)
self.frame_sequence_cache.add("frame_sequence_potential_energy", [], single=False, update=True)
self.frame_sequence_cache.add("frame_sequence_kinetic_energy", [], single=False, update=True)
self.frame_sequence_cache.add("frame_sequence_temperature", [], single=False, update=True)
self.frame_sequence_cache.add("frame_sequence_time", [], single=False, update=True)
self.frame_sequence_cache.add("frame_sequence_to_sampling_ref", single=False, update=True)
# Cache for storing system information
self.system_cache = CacheService(self.parser_context)
......@@ -110,10 +112,9 @@ class NWChemMainParser(MainHierarchicalParser):
def geo_opt_module(self):
return SM( " NWChem Geometry Optimization",
sections=["section_method", "section_frame_sequence", "section_sampling_method", "x_nwchem_section_geo_opt_module"],
onClose={
"section_sampling_method": self.save_geo_opt_sampling_id,
"section_frame_sequence": self.save_local_frames_ref,
},
# onClose={
# "section_sampling_method": self.save_geo_opt_sampling_id,
# },
subFlags=SM.SubFlags.Sequenced,
subMatchers=[
SM( r" maximum gradient threshold \(gmax\) =\s+(?P<geometry_optimization_threshold_force__forceAu>{})".format(self.regexs.float)),
......@@ -197,12 +198,12 @@ class NWChemMainParser(MainHierarchicalParser):
SM(" DFT ENERGY GRADIENTS"),
SM(" QMD Run Information",
subMatchers=[
SM(" Time elapsed \(fs\) :(?P<x_nwchem_qmd_step_time__fs>{})".format(self.regexs.float)),
SM(" Time elapsed \(fs\) :\s+(?P<x_nwchem_qmd_step_time__fs>{})".format(self.regexs.float)),
SM(" Kin. energy \(a\.u\.\):\s+{}\s+(?P<x_nwchem_qmd_step_kinetic_energy__hartree>{})".format(self.regexs.int, self.regexs.float)),
SM(" Pot. energy \(a\.u\.\):\s+{}\s+(?P<x_nwchem_qmd_step_potential_energy__hartree>{})".format(self.regexs.int, self.regexs.float)),
SM(" Tot. energy \(a\.u\.\):\s+{}\s+(?P<x_nwchem_qmd_step_total_energy__hartree>{})".format(self.regexs.int, self.regexs.float)),
SM(" Target temp. (K) :\s+{}\s+(?P<x_nwchem_qmd_step_target_temperature__K>{})".format(self.regexs.int, self.regexs.float)),
SM(" Current temp. (K) :\s+{}\s+(?P<x_nwchem_qmd_step_temperature__K>{})".format(self.regexs.int, self.regexs.float)),
SM(" Target temp\. \(K\) :\s+{}\s+(?P<x_nwchem_qmd_step_target_temperature__K>{})".format(self.regexs.int, self.regexs.float)),
SM(" Current temp\. \(K\) :\s+{}\s+(?P<x_nwchem_qmd_step_temperature__K>{})".format(self.regexs.int, self.regexs.float)),
SM(" Dipole \(a\.u\.\) :\s+{0}\s+({1}\s+{1}\s+{1})".format(self.regexs.int, self.regexs.float), startReTransform=self.dipole_transform)
]
)
......@@ -347,12 +348,26 @@ class NWChemMainParser(MainHierarchicalParser):
def onClose_section_frame_sequence(self, backend, gIndex, section):
self.frame_sequence_cache.addValue("number_of_frames_in_sequence")
self.frame_sequence_cache.addArrayValues("frame_sequence_local_frames_ref")
self.frame_sequence_cache.addValue("frame_sequence_to_sampling_ref")
potential_energy = np.array(self.frame_sequence_cache["frame_sequence_potential_energy"])
backend.addArrayValues("frame_sequence_potential_energy", potential_energy)
if potential_energy.size != 0:
backend.addArrayValues("frame_sequence_potential_energy", potential_energy)
backend.addArrayValues("frame_sequence_potential_energy_stats", np.array([potential_energy.mean(), potential_energy.std()]))
kin_energy = np.array(self.frame_sequence_cache["frame_sequence_kinetic_energy"])
backend.addArrayValues("frame_sequence_kinetic_energy", kin_energy)
if kin_energy.size != 0:
backend.addArrayValues("frame_sequence_kinetic_energy", kin_energy)
backend.addArrayValues("frame_sequence_kinetic_energy_stats", np.array([kin_energy.mean(), kin_energy.std()]))
temp = np.array(self.frame_sequence_cache["frame_sequence_temperature"])
if temp.size != 0:
backend.addArrayValues("frame_sequence_temperature", temp)
backend.addArrayValues("frame_sequence_temperature_stats", np.array([temp.mean(), temp.std()]))
time = np.array(self.frame_sequence_cache["frame_sequence_time"])
if time.size != 0:
backend.addArrayValues("frame_sequence_time", time)
self.frame_sequence_cache.clear()
......@@ -373,6 +388,12 @@ class NWChemMainParser(MainHierarchicalParser):
kin_energy = section.get_latest_value("x_nwchem_qmd_step_kinetic_energy")
self.frame_sequence_cache["frame_sequence_kinetic_energy"].append(kin_energy)
temp = section.get_latest_value("x_nwchem_qmd_step_temperature")
self.frame_sequence_cache["frame_sequence_temperature"].append(temp)
time = section.get_latest_value("x_nwchem_qmd_step_time")
self.frame_sequence_cache["frame_sequence_time"].append(time)
def onClose_x_nwchem_section_geo_opt_step(self, backend, gIndex, section):
self.frame_sequence_cache["number_of_frames_in_sequence"] += 1
pot_ener = section.get_latest_value("x_nwchem_geo_opt_step_energy")
......@@ -384,6 +405,9 @@ class NWChemMainParser(MainHierarchicalParser):
self.method_index = gIndex
self.save_method = True
def onOpen_section_sampling_method(self, backend, gIndex, section):
self.frame_sequence_cache["frame_sequence_to_sampling_ref"] = gIndex
def onOpen_x_nwchem_section_qmd_module(self, backend, gIndex, section):
self.sampling_method_cache["sampling_method"] = "molecular_dynamics"
......@@ -596,10 +620,6 @@ class NWChemMainParser(MainHierarchicalParser):
def add_frame_reference(self, backend, gIndex, section):
self.frame_sequence_cache["frame_sequence_local_frames_ref"].append(gIndex)
def save_local_frames_ref(self, backend, gIndex, section):
backend.addArrayValues("frame_sequence_local_frames_ref", np.array(self.frame_sequence_local_frames_ref))
self.frame_sequence_local_frames_ref = []
#=======================================================================
# Start match transforms
def dipole_transform(self, backend, groups):
......
......@@ -16,4 +16,8 @@ object NWChemParserSpec extends Specification {
"test geo_opt with json" >> {
ParserRun.parse(NWChemParser, "parsers/nwchem/test/examples/geo_opt/output.out", "json") must_== ParseResult.ParseSuccess
}
"test md with json" >> {
ParserRun.parse(NWChemParser, "parsers/nwchem/test/examples/md/output.out", "json") must_== ParseResult.ParseSuccess
}
}
start qmd_dft_h2o_svr
echo
print low
geometry noautosym noautoz
O 0.00000000 -0.01681748 0.11334792
H 0.00000000 0.81325914 -0.34310308
H 0.00000000 -0.67863597 -0.56441201
end
basis
* library 6-31G*
end
dft
xc pbe0
end
qmd
nstep_nucl 5
dt_nucl 10.0
targ_temp 200.0
com_step 10
thermostat svr 100.0
print_xyz 5
end
task dft qmd
argument 1 = input.nw
============================== echo of input deck ==============================
start qmd_dft_h2o_svr
echo
print low
geometry noautosym noautoz
O 0.00000000 -0.01681748 0.11334792
H 0.00000000 0.81325914 -0.34310308
H 0.00000000 -0.67863597 -0.56441201
end
basis
* library 6-31G*
end
dft
xc pbe0
end
qmd
nstep_nucl 5
dt_nucl 10.0
targ_temp 200.0
com_step 10
thermostat svr 100.0
print_xyz 5
end
task dft qmd
================================================================================
Northwest Computational Chemistry Package (NWChem) 6.6
------------------------------------------------------
Environmental Molecular Sciences Laboratory
Pacific Northwest National Laboratory
Richland, WA 99352
Copyright (c) 1994-2015
Pacific Northwest National Laboratory
Battelle Memorial Institute
NWChem is an open-source computational chemistry package
distributed under the terms of the
Educational Community License (ECL) 2.0
A copy of the license is included with this distribution
in the LICENSE.TXT file
ACKNOWLEDGMENT
--------------
This software and its documentation were developed at the
EMSL at Pacific Northwest National Laboratory, a multiprogram
national laboratory, operated for the U.S. Department of Energy
by Battelle under Contract Number DE-AC05-76RL01830. Support
for this work was provided by the Department of Energy Office
of Biological and Environmental Research, Office of Basic
Energy Sciences, and the Office of Advanced Scientific Computing.
Job information
---------------
hostname = lenovo700
program = nwchem
date = Wed Sep 14 15:04:49 2016
compiled = Mon_Feb_15_08:24:17_2016
source = /build/nwchem-MF0R1k/nwchem-6.6+r27746
nwchem branch = 6.6
nwchem revision = 27746
ga revision = 10594
input = input.nw
prefix = qmd_dft_h2o_svr.
data base = ./qmd_dft_h2o_svr.db
status = startup
nproc = 1
time left = -1s
Memory information
------------------
heap = 13107198 doubles = 100.0 Mbytes
stack = 13107195 doubles = 100.0 Mbytes
global = 26214400 doubles = 200.0 Mbytes (distinct from heap & stack)
total = 52428793 doubles = 400.0 Mbytes
verify = yes
hardfail = no
Directory information
---------------------
0 permanent = .
0 scratch = .
NWChem Input Module
-------------------
Scaling coordinates for geometry "geometry" by 1.889725989
(inverse scale = 0.529177249)
Geometry "geometry" -> ""
-------------------------
Output coordinates in angstroms (scale by 1.889725989 to convert to a.u.)
No. Tag Charge X Y Z
---- ---------------- ---------- -------------- -------------- --------------
1 O 8.0000 0.00000000 -0.01682581 0.11342109
2 H 1.0000 0.00000000 0.81325081 -0.34302991
3 H 1.0000 0.00000000 -0.67864430 -0.56433884
Atomic Mass
-----------
O 15.994910
H 1.007825
Effective nuclear repulsion energy (a.u.) 9.2887672039
Nuclear Dipole moment (a.u.)
----------------------------
X Y Z
---------------- ---------------- ----------------
0.0000000000 0.0000000000 0.0000000000
XYZ format geometry
-------------------
3
geometry
O 0.00000000 -0.01682581 0.11342109
H 0.00000000 0.81325081 -0.34302991
H 0.00000000 -0.67864430 -0.56433884
==============================================================================
internuclear distances
------------------------------------------------------------------------------
center one | center two | atomic units | angstroms
------------------------------------------------------------------------------
2 H | 1 O | 1.79013 | 0.94730
3 H | 1 O | 1.79012 | 0.94729
------------------------------------------------------------------------------
number of included internuclear distances: 2
==============================================================================
==============================================================================
internuclear angles
------------------------------------------------------------------------------
center 1 | center 2 | center 3 | degrees
------------------------------------------------------------------------------
2 H | 1 O | 3 H | 105.51
------------------------------------------------------------------------------
number of included internuclear angles: 1
==============================================================================
library name resolved from: .nwchemrc
library file name is: </home/lauri/nwchem-6.6/src/basis/libraries/>
Summary of "ao basis" -> "" (cartesian)
------------------------------------------------------------------------------
Tag Description Shells Functions and Types
---------------- ------------------------------ ------ ---------------------
* 6-31G* on all atoms
calling qmd_driver
NWChem QMD Module
-----------------
QMD Run Parameters
------------------
No. of nuclear steps: 5
Nuclear time step: 10.00
Target temp. (K): 200.00
Thermostat: svr
Tau: 100.00
Random seed: -8677
Nuclear integrator: velocity-verlet
No restart file found
Beginning with random velocities
Current temp. (K): 830.59
NWChem DFT Module
-----------------
Basis "ao basis" -> "ao basis" (cartesian)
-----
O (Oxygen)
----------
Exponent Coefficients
-------------- ---------------------------------------------------------
1 S 5.48467170E+03 0.001831
1 S 8.25234950E+02 0.013950
1 S 1.88046960E+02 0.068445
1 S 5.29645000E+01 0.232714
1 S 1.68975700E+01 0.470193
1 S 5.79963530E+00 0.358521
2 S 1.55396160E+01 -0.110778
2 S 3.59993360E+00 -0.148026
2 S 1.01376180E+00 1.130767
3 P 1.55396160E+01 0.070874
3 P 3.59993360E+00 0.339753
3 P 1.01376180E+00 0.727159
4 S 2.70005800E-01 1.000000
5 P 2.70005800E-01 1.000000
6 D 8.00000000E-01 1.000000
H (Hydrogen)
------------
Exponent Coefficients
-------------- ---------------------------------------------------------
1 S 1.87311370E+01 0.033495
1 S 2.82539370E+00 0.234727
1 S 6.40121700E-01 0.813757
2 S 1.61277800E-01 1.000000
Summary of "ao basis" -> "ao basis" (cartesian)
------------------------------------------------------------------------------
Tag Description Shells Functions and Types
---------------- ------------------------------ ------ ---------------------
O 6-31G* 6 15 3s2p1d
H 6-31G* 2 2 2s
Caching 1-el integrals
Time after variat. SCF: 0.1
Time prior to 1st pass: 0.1
Total DFT energy = -76.325066149291
One electron energy = -123.272247284013
Coulomb energy = 46.936105040748
Exchange-Corr. energy = -9.277691109951
Nuclear repulsion energy = 9.288767203925
Numeric. integr. density = 10.000001227433
Total iterative time = 0.4s
DFT ENERGY GRADIENTS
atom coordinates gradient
x y z x y z
1 O 0.000000 -0.031796 0.214335 0.000000 0.003514 -0.023728
2 H 0.000000 1.536821 -0.648233 -0.000000 -0.014425 0.009983
3 H 0.000000 -1.282452 -1.066446 0.000000 0.010912 0.013746
NWChem DFT Module
-----------------
Caching 1-el integrals
Time after variat. SCF: 0.6
Time prior to 1st pass: 0.6
Total DFT energy = -76.324876801053
One electron energy = -123.293471622871
Coulomb energy = 46.945910580941
Exchange-Corr. energy = -9.278925089848
Nuclear repulsion energy = 9.301609330725
Numeric. integr. density = 10.000001001327
Total iterative time = 0.4s
DFT ENERGY GRADIENTS
atom coordinates gradient
x y z x y z
1 O -0.000000 -0.030541 0.214843 0.000000 0.017462 -0.022827
2 H 0.000000 1.520797 -0.645137 -0.000000 -0.023895 0.014570
3 H -0.000000 -1.286347 -1.077605 0.000000 0.006434 0.008257
QMD Run Information
-------------------
Time elapsed (fs) : 0.241888
Kin. energy (a.u.): 1 0.003270
Pot. energy (a.u.): 1 -76.324877
Tot. energy (a.u.): 1 -76.321607
Target temp. (K) : 1 200.00
Current temp. (K) : 1 688.40
Dipole (a.u.) : 1 -1.297680E-11 1.171852E-01 -8.161034E-01
NWChem DFT Module
-----------------
Caching 1-el integrals
Time after variat. SCF: 1.1
Time prior to 1st pass: 1.1
Total DFT energy = -76.324529295229
One electron energy = -123.304435531451
Coulomb energy = 46.950661720540
Exchange-Corr. energy = -9.279572169719
Nuclear repulsion energy = 9.308816685402
Numeric. integr. density = 10.000000736576
Total iterative time = 0.4s
DFT ENERGY GRADIENTS
atom coordinates gradient
x y z x y z
1 O -0.000000 -0.029445 0.215385 0.000000 0.029824 -0.021815
2 H 0.000000 1.507317 -0.643053 -0.000000 -0.032130 0.018615
3 H -0.000000 -1.290264 -1.088292 0.000000 0.002306 0.003199
QMD Run Information
-------------------
Time elapsed (fs) : 0.483777
Kin. energy (a.u.): 2 0.003055
Pot. energy (a.u.): 2 -76.324529
Tot. energy (a.u.): 2 -76.321474
Target temp. (K) : 2 200.00
Current temp. (K) : 2 643.13
Dipole (a.u.) : 2 -2.348946E-11 1.141435E-01 -8.184600E-01
NWChem DFT Module
-----------------