Commit 52e6a24d authored by Lauri Himanen's avatar Lauri Himanen
Browse files

The optimised geometry is now read from an external file with ASE.

parent 48047674
......@@ -4,7 +4,6 @@ import logging
from nomadcore.simple_parser import SimpleMatcher as SM
from nomadcore.simple_parser import extractOnCloseTriggers
from nomadcore.caching_backend import CachingLevel
from nomadcore.baseclasses import CacheMode
from inputparser import CP2KInputParser
logger = logging.getLogger("nomad")
......@@ -48,11 +47,11 @@ class CommonMatcher(object):
#=======================================================================
# Cached values
self.cache_service.add_cache_object("simulation_cell", CacheMode.SINGLE_IN_MULTI_OUT)
self.cache_service.add_cache_object("number_of_scf_iterations", CacheMode.MULTI_IN_MULTI_OUT, 0)
self.cache_service.add_cache_object("atom_positions", CacheMode.MULTI_IN_MULTI_OUT)
self.cache_service.add_cache_object("atom_labels", CacheMode.SINGLE_IN_MULTI_OUT)
self.cache_service.add_cache_object("number_of_atoms", CacheMode.SINGLE_IN_MULTI_OUT)
self.cache_service.add_cache_object("simulation_cell", single=False, update=False)
self.cache_service.add_cache_object("number_of_scf_iterations", 0)
self.cache_service.add_cache_object("atom_positions", single=False, update=True)
self.cache_service.add_cache_object("atom_labels", single=False, update=False)
self.cache_service.add_cache_object("number_of_atoms", single=False, update=False)
#===========================================================================
# SimpleMatchers
......@@ -259,9 +258,6 @@ class CommonMatcher(object):
def onClose_section_single_configuration_calculation(self, backend, gIndex, section):
"""
"""
self.cache_service.push_value("number_of_scf_iterations")
self.cache_service["number_of_scf_iterations"] = 0
# Write the references to section_method and section_system
backend.addValue('single_configuration_to_calculation_method_ref', self.section_method_index)
backend.addValue('single_configuration_calculation_to_system_ref', self.section_system_index)
......
from nomadcore.simple_parser import SimpleMatcher as SM
from nomadcore.baseclasses import MainHierarchicalParser, CacheMode
from nomadcore.baseclasses import MainHierarchicalParser
from commonmatcher import CommonMatcher
from nomadcore.caching_backend import CachingLevel
import logging
import ase.io
logger = logging.getLogger("nomad")
......@@ -16,11 +17,12 @@ class CP2KGeoOptParser(MainHierarchicalParser):
"""
super(CP2KGeoOptParser, self).__init__(file_path, parser_context)
self.setup_common_matcher(CommonMatcher(parser_context))
self.traj_iterator = None
#=======================================================================
# Cached values
self.cache_service.add_cache_object("number_of_frames_in_sequence", CacheMode.MULTI_IN_MULTI_OUT, 0)
self.cache_service.add_cache_object("frame_sequence_potential_energy", CacheMode.MULTI_IN_MULTI_OUT, [])
self.cache_service.add_cache_object("number_of_frames_in_sequence", 0, single=True, update=True)
self.cache_service.add_cache_object("frame_sequence_potential_energy", [], single=True, update=True)
#=======================================================================
# Cache levels
......@@ -34,7 +36,8 @@ class CP2KGeoOptParser(MainHierarchicalParser):
" *** STARTING GEOMETRY OPTIMIZATION ***".replace("*", "\*"),
sections=["section_frame_sequence", "section_sampling_method"],
subMatchers=[
SM( " REQUESTED STRUCTURE DATA",
SM( " -------- Informations at step =\s+{}\s+------------".format(self.cm.regex_i),
forwardMatch=True,
name="geooptstep",
repeats=True,
sections=["section_single_configuration_calculation", "section_system"],
......@@ -78,7 +81,7 @@ class CP2KGeoOptParser(MainHierarchicalParser):
subMatchers=[
self.cm.header(),
self.cm.quickstep(),
]
],
),
self.geo_opt
]
......@@ -94,6 +97,15 @@ class CP2KGeoOptParser(MainHierarchicalParser):
energy = section["x_cp2k_optimization_energy"][0]
self.cache_service["frame_sequence_potential_energy"].append(energy)
def onClose_section_method(self, backend, gIndex, section):
traj_file = self.file_service.get_file_by_id("trajectory")
try:
if traj_file is not None:
self.traj_iterator = ase.io.iread(traj_file)
except ValueError:
# The format was not supported by ase
pass
#===========================================================================
# adHoc functions
def adHoc_geo_opt_converged(self):
......@@ -111,8 +123,26 @@ class CP2KGeoOptParser(MainHierarchicalParser):
return wrapper
def adHoc_step(self):
"""Called when the geometry optimization did not converge.
"""Called when all the step information has been retrieved from the
output file. Here further information is gathered from external files.
"""
def wrapper(parser):
self.cache_service["number_of_frames_in_sequence"] += 1
# Get the next position from the trajectory file
if self.traj_iterator is not None:
atoms = next(self.traj_iterator)
pos = atoms.positions
self.cache_service["atom_positions"] = pos
return wrapper
def adHoc_setup_traj_file(self):
def wrapper(parser):
print "HERE"
return wrapper
def debug(self):
def wrapper(parser):
print "FOUND"
return wrapper
......@@ -3,7 +3,7 @@ import re
import logging
import cPickle as pickle
import numpy as np
from nomadcore.baseclasses import BasicParser, CacheMode
from nomadcore.baseclasses import BasicParser
from cp2kparser.generic.inputparsing import *
logger = logging.getLogger("nomad")
......@@ -40,9 +40,11 @@ class CP2KInputParser(BasicParser):
self.input_tree = None
self.input_lines = None
self.force_file_name = None
self.trajectory_file_name = ""
# Declare the cached values here
self.cache_service.add_cache_object("configuration_periodic_dimensions", CacheMode.SINGLE_IN_MULTI_OUT)
#=======================================================================
# Cached values
self.cache_service.add_cache_object("configuration_periodic_dimensions", single=False, update=False)
def parse(self):
......@@ -143,11 +145,11 @@ class CP2KInputParser(BasicParser):
#=======================================================================
# Single point force file name
# force_file = self.input_tree.get_keyword("FORCE_EVAL/PRINT/FORCES/FILENAME")
force_file = self.force_file_name
if force_file is not None and force_file != "__STD_OUT__":
force_file_path = self.normalize_x_cp2k_path(force_file, "xyz")
self.file_service.set_file_id(force_file_path, "force_file_single_point")
self.setup_force_file_name()
#=======================================================================
# Trajectory file name
self.setup_trajectory_file_name()
#=======================================================================
# Stress tensor calculation method
......@@ -163,21 +165,40 @@ class CP2KInputParser(BasicParser):
if stress_tensor_method is not None:
self.backend.addValue("stress_tensor_method", stress_tensor_method)
def normalize_x_cp2k_path(self, path, extension, name=""):
def normalize_x_cp2k_path(self, path):
"""The paths in CP2K input can be given in many ways. This function
tries to normalize these forms into a valid path.
"""
if name:
name = "-" + name
project_name = self.input_tree.get_keyword("GLOBAL/PROJECT_NAME")
# Path is exactly as given
if path.startswith("="):
normalized_path = path[1:]
# Path is relative, no project name added
elif re.match(r"./", path):
normalized_path = "{}{}-1_0.{}".format(path, name, extension)
normalized_path = path
# Path is relative, project name added
else:
normalized_path = "{}-{}{}-1_0.{}".format(project_name, path, name, extension)
project_name = self.input_tree.get_keyword("GLOBAL/PROJECT_NAME")
normalized_path = "{}-{}".format(project_name, path)
return normalized_path
def setup_force_file_name(self):
"""Setup the force file path.
"""
force_file = self.force_file_name
extension = "xyz"
if force_file is not None and force_file != "__STD_OUT__":
normalized_path = self.normalize_x_cp2k_path(self.force_file_name)
final_path = "{}-1_0.{}".format(normalized_path, extension)
self.file_service.set_file_id(final_path, "force_file_single_point")
def setup_trajectory_file_name(self):
"""Setup the trajectory file path.
"""
extension = "xyz"
normalized_path = self.normalize_x_cp2k_path(self.trajectory_file_name)
final_path = "{}pos-1.{}".format(normalized_path, extension)
self.file_service.set_file_id(final_path, "trajectory")
def fill_input_tree(self, file_path):
"""Parses a CP2K input file into an object tree.
......@@ -256,6 +277,9 @@ class CP2KInputParser(BasicParser):
if path == "FORCE_EVAL/PRINT/FORCES":
if keyword_name == "FILENAME":
self.force_file_name = keyword_value
if path == "MOTION/PRINT/TRAJECTORY":
if keyword_name == "FILENAME":
self.trajectory_file_name = keyword_value
def fill_metadata(self):
"""Goes through the input data and pushes everything to the
......
......@@ -47,5 +47,11 @@ class CP2KSinglePointParser(MainHierarchicalParser):
else:
logger.warning("The file containing the forces printed by ENERGY_FORCE calculation could not be found.")
# Only in the single configuration calculations the number of scf
# iterations is given. E.g. in geometry optimization there are multiple
# scf calculations so this loses it's meaning sort of.
self.cache_service.push_value("number_of_scf_iterations")
self.cache_service["number_of_scf_iterations"] = 0
#===========================================================================
# adHoc functions
NONBONDED NEIGHBOR LISTS IN angstrom (PROCESS 0)
Atom-A X Y Z Atom-B X Y Z Cell(i,j,k) Distance ONFO VDW-scale EI-scale
3 -0.491324 1.573799 -2.426806 1 -0.178478 1.376642 -1.543920 0 0 0 0.9572
1 -0.178478 1.376642 -1.543920 2 0.001339 2.233125 -1.156189 0 0 0 0.9572
3 -0.491324 1.573799 -2.426806 2 0.001339 2.233125 -1.156189 0 0 0 1.5139
Total number of neighbor interactions for process 0: 3
3
i = 1, E = -17.1643447508
O 12.2353220000 1.3766420000 10.8698800000
H 12.4189404412 2.2491702452 11.2652097702
H 11.9149490419 1.5748472981 9.9693615826
3
i = 2, E = -17.1643578234
O 12.2353220000 1.3766420000 10.8698800000
H 12.4200101767 2.2468980462 11.2643381000
H 11.9115895109 1.5703371433 9.9671556911
3
i = 3, E = -17.1643766833
O 12.2353220000 1.3766420000 10.8698800000
H 12.4233662348 2.2477209582 11.2631307672
H 11.9069167231 1.5604096790 9.9689943947
3
i = 4, E = -17.1644565615
O 12.2353220000 1.3766420000 10.8698800000
H 12.4423251455 2.2644837735 11.2309777118
H 11.8922856584 1.5514786478 9.9754905025
3
i = 5, E = -17.1645203350
O 12.2353220000 1.3766420000 10.8698800000
H 12.4431348352 2.2577511412 11.2267313861
H 11.8919233389 1.5518052357 9.9721549140
3
i = 6, E = -17.1646190312
O 12.2353220000 1.3766420000 10.8698800000
H 12.4736086753 2.2567157451 11.2132013565
H 11.9047281262 1.5469224568 9.9655629396
3
i = 7, E = -17.1646198446
O 12.2353220000 1.3766420000 10.8698800000
H 12.4742379375 2.2563137117 11.2125606428
H 11.9051180271 1.5469372555 9.9662564829
3
i = 8, E = -17.1646203294
O 12.2353220000 1.3766420000 10.8698800000
H 12.4749724176 2.2565225180 11.2119388324
H 11.9039582865 1.5481084051 9.9667841313
3
i = 9, E = -17.1646204547
O 12.2353220000 1.3766420000 10.8698800000
H 12.4752917518 2.2561254118 11.2121163217
H 11.9031806595 1.5481019598 9.9671120050
3
i = 10, E = -17.1646204237
O 12.2353220000 1.3766420000 10.8698800000
H 12.4755864059 2.2562726297 11.2125118440
H 11.9026358443 1.5480115451 9.9673498890
3
i = 11, E = -17.1646204766
O 12.2353220000 1.3766420000 10.8698800000
H 12.4754906682 2.2560912723 11.2123985687
H 11.9026556670 1.5480524898 9.9673305122
3
i = 12, E = -17.1646347711
O 12.2353220000 1.3766420000 10.8698800000
H 12.4754916182 2.2560930719 11.2123996927
H 11.9026554703 1.5480520836 9.9673307045
&GLOBAL
PROJECT H2O
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&SUBSYS
&CELL
ABC 12.4138 12.4138 12.4138
&END CELL
&COORD
O 12.235322 1.376642 10.869880
H 12.415139 2.233125 11.257611
H 11.922476 1.573799 9.986994
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME ../BASIS_SET
POTENTIAL_FILE_NAME ../GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-7
&END QS
&MGRID
CUTOFF 200
NGRIDS 4
REL_CUTOFF 30
&END MGRID
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-05
MAX_SCF 200
&DIAGONALIZATION T
ALGORITHM STANDARD
&END DIAGONALIZATION
&MIXING T
ALPHA 0.5
METHOD PULAY_MIXING
NPULAY 5
&END MIXING
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&END FORCE_EVAL
&MOTION
&GEO_OPT
TYPE MINIMIZATION
MAX_DR 1.0E-03
MAX_FORCE 1.0E-03
RMS_DR 1.0E-03
RMS_FORCE 1.0E-03
MAX_ITER 200
OPTIMIZER CG
&CG
MAX_STEEP_STEPS 0
RESTART_LIMIT 9.0E-01
&END CG
&END GEO_OPT
&CONSTRAINT
&FIXED_ATOMS
COMPONENTS_TO_FIX XYZ
LIST 1
&END FIXED_ATOMS
&END CONSTRAINT
&END MOTION
This diff is collapsed.
......@@ -438,7 +438,7 @@ class TestPreprocessor(unittest.TestCase):
result = get_result("input_preprocessing/variable", "x_cp2k_CP2K_INPUT.GLOBAL.PROJECT_NAME", optimize=False)
self.assertEqual(result, "variable_test")
def test_variable_mutiple(self):
def test_variable_multiple(self):
result = get_result("input_preprocessing/variable_multiple", "x_cp2k_CP2K_INPUT.FORCE_EVAL.DFT.MGRID.CUTOFF", optimize=False)
self.assertEqual(result, 50)
......@@ -448,15 +448,39 @@ class TestGeoOpt(unittest.TestCase):
@classmethod
def setUpClass(cls):
cls.results = get_results("geo_opt", "section_run")
cls.results = get_results("geo_opt/cg", "section_run")
def test_geometry_optimization_converged(self):
result = self.results["geometry_optimization_converged"]
self.assertTrue(result)
def test_number(self):
def test_number_of_frames_in_sequence(self):
result = self.results["number_of_frames_in_sequence"]
self.assertEqual(result, 12)
self.assertEqual(result, 7)
def test_atom_positions(self):
result = self.results["atom_positions"]
expected_start = convert_unit(
np.array([
[12.2353220000, 1.3766420000, 10.8698800000],
[12.4175775999, 2.2362362573, 11.2616216864],
[11.9271436933, 1.5723516602, 10.0115134757],
]),
"angstrom"
)
expected_end = convert_unit(
np.array([
[12.2353220000, 1.3766420000, 10.8698800000],
[12.4958164689, 2.2307248873, 11.3354322515],
[11.9975558616, 1.5748085240, 10.0062792262],
]),
"angstrom"
)
result_start = result[0,:,:]
result_end = result[-1,:,:]
self.assertTrue(np.array_equal(result_start, expected_start))
self.assertTrue(np.array_equal(result_end, expected_end))
#===============================================================================
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment