Skip to content

Commit

Permalink
Stop using PDB because its precision is too low
Browse files Browse the repository at this point in the history
  • Loading branch information
oschuett committed Nov 24, 2017
1 parent be2989c commit 96f06d2
Show file tree
Hide file tree
Showing 6 changed files with 166 additions and 37 deletions.
7 changes: 3 additions & 4 deletions aiida_cp2k/calculations/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ def _init_internal_params(self):
self._PROJECT_NAME = 'aiida'
self._RESTART_FILE_NAME = self._PROJECT_NAME + '-1.restart'
self._PARENT_CALC_FOLDER_NAME = 'parent_calc/'
self._COORDS_FILE_NAME = 'aiida.coords.pdb'
self._COORDS_FILE_NAME = 'aiida.coords.xyz'
self._default_parser = 'cp2k'

# --------------------------------------------------------------------------
Expand Down Expand Up @@ -109,16 +109,15 @@ def _prepare_for_submission(self, tempfolder, inputdict):
# write cp2k input file
inp = Cp2kInput(params)
inp.add_keyword("GLOBAL/PROJECT", self._PROJECT_NAME)
inp.add_keyword("MOTION/PRINT/TRAJECTORY/FORMAT", "PDB")
if structure is not None:
struct_fn = tempfolder.get_abs_path(self._COORDS_FILE_NAME)
structure.get_ase().write(struct_fn, format="proteindatabank")
structure.export(struct_fn, fileformat="xyz")
for i, a in enumerate('ABC'):
val = '{:<15} {:<15} {:<15}'.format(*structure.cell[i])
inp.add_keyword('FORCE_EVAL/SUBSYS/CELL/'+a, val)
topo = "FORCE_EVAL/SUBSYS/TOPOLOGY"
inp.add_keyword(topo + "/COORD_FILE_NAME", self._COORDS_FILE_NAME)
inp.add_keyword(topo + "/COORD_FILE_FORMAT", "pdb")
inp.add_keyword(topo + "/COORD_FILE_FORMAT", "XYZ")
inp_fn = tempfolder.get_abs_path(self._INPUT_FILE_NAME)
with open(inp_fn, "w") as f:
f.write(inp.render())
Expand Down
1 change: 1 addition & 0 deletions test/run_tests.sh
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ verdi daemon start
./test_no_struct.py cp2k@localhost
./test_restart.py cp2k@localhost
./test_failure.py cp2k@localhost
./test_precision.py cp2k@localhost

echo "All tests have passed :-)"
#EOF
2 changes: 1 addition & 1 deletion test/test_dft.py
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@
wait_for_calc(calc)

# check energy
expected_energy = -17.1566368539
expected_energy = -17.1566361119
if abs(calc.res.energy - expected_energy) < 1e-10:
print("OK, energy has the expected value")
else:
Expand Down
4 changes: 2 additions & 2 deletions test/test_geopt.py
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@
assert calc.res.exceeded_walltime is False

# check energy
expected_energy = -1.14009973333
expected_energy = -1.14009973178
if abs(calc.res.energy - expected_energy) < 1e-10:
print("OK, energy has the expected value")
else:
Expand All @@ -110,7 +110,7 @@
sys.exit(3)

# check geometry
expected_dist = 0.736125210983
expected_dist = 0.736103879818
dist = calc.out.output_structure.get_ase().get_distance(0, 1)
if abs(dist - expected_dist) < 1e-7:
print("OK, H-H distance has the expected value")
Expand Down
71 changes: 41 additions & 30 deletions test/test_mm.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,37 @@
# calc object
calc = code.new_calc()

# force field
with open("/tmp/water.pot", "w") as f:
f.write("""BONDS
H H 0.000 1.5139
O H 450.000 0.9572
ANGLES
H O H 55.000 104.5200
DIHEDRALS
IMPROPER
NONBONDED
H 0.000000 -0.046000 0.224500
O 0.000000 -0.152100 1.768200
HBOND CUTHB 0.5
END""")
water_pot = SinglefileData(file="/tmp/water.pot")
calc.use_file(water_pot, linkname="water_pot")

# structure using pdb format, because it also carries topology information
atoms = ase.build.molecule('H2O')
atoms.center(vacuum=10.0)
atoms.write("/tmp/coords.pdb", format="proteindatabank")
cell = atoms.cell
coords_pdb = SinglefileData(file="/tmp/coords.pdb")
calc.use_file(coords_pdb, linkname="coords_pdb")

# parameters
# based on cp2k/tests/Fist/regtest-1-1/water_1.inp
parameters = ParameterData(dict={
Expand All @@ -57,7 +88,16 @@
'GMAX': 24,
'O_SPLINE': 6
}}
}
},
'SUBSYS': {
'CELL': {
'ABC': '%f %f %f' % tuple(atoms.cell.diagonal()),
},
'TOPOLOGY': {
'COORD_FILE_NAME': 'coords.pdb',
'COORD_FILE_FORMAT': 'PDB',
},
},
},
'GLOBAL': {
'CALLGRAPH': 'master',
Expand All @@ -66,35 +106,6 @@
})
calc.use_parameters(parameters)

# force field
with open("/tmp/water.pot", "w") as f:
f.write("""BONDS
H H 0.000 1.5139
O H 450.000 0.9572
ANGLES
H O H 55.000 104.5200
DIHEDRALS
IMPROPER
NONBONDED
H 0.000000 -0.046000 0.224500
O 0.000000 -0.152100 1.768200
HBOND CUTHB 0.5
END""")
water_pot = SinglefileData(file="/tmp/water.pot")
calc.use_file(water_pot, linkname="water_pot")

# structure
atoms = ase.build.molecule('H2O')
atoms.center(vacuum=10.0)
structure = StructureData(ase=atoms)
calc.use_structure(structure)

# settings
settings_dict = {'additional_retrieve_list': ["runtime.callgraph"]}
settings = ParameterData(dict=settings_dict)
Expand Down
118 changes: 118 additions & 0 deletions test/test_precision.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,118 @@
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
###########################################################################
# Copyright (c), The AiiDA team. All rights reserved. #
# This file is part of the AiiDA code. #
# #
# The code is hosted on GitHub at https://github.com/cp2k/aiida-cp2k #
# For further information on the license, see the LICENSE.txt file #
# For further information please visit http://www.aiida.net #
###########################################################################
from __future__ import print_function

import sys
import ase.build
import numpy as np
from utils import wait_for_calc

from aiida import load_dbenv, is_dbenv_loaded
from aiida.backends import settings
if not is_dbenv_loaded():
load_dbenv(profile=settings.AIIDADB_PROFILE)

from aiida.common.example_helpers import test_and_get_code # noqa
from aiida.orm.data.structure import StructureData # noqa
from aiida.orm.data.parameter import ParameterData # noqa
from aiida.orm.data.singlefile import SinglefileData # noqa


# ==============================================================================
if len(sys.argv) != 2:
print("Usage: test_precision.py <code_name>")
sys.exit(1)

codename = sys.argv[1]
code = test_and_get_code(codename, expected_code_type='cp2k')

print("Testing structure roundtrip precision ase->aiida->cp2k->aiida->ase...")

# calc object
calc = code.new_calc()

# structure
epsilon = 1e-10 # expected precision in Angstrom
dist = 0.74 + epsilon
atoms = ase.Atoms('H2', positions=[(0, 0, 0), (0, 0, dist)], cell=[4, 4, 4])
structure = StructureData(ase=atoms)
calc.use_structure(structure)


# parameters
parameters = ParameterData(dict={
'GLOBAL': {
'RUN_TYPE': 'MD',
},
'MOTION': {
'MD': {
'TIMESTEP': 0.0, # do not move atoms
'STEPS': 1,
},
},
'FORCE_EVAL': {
'METHOD': 'Quickstep',
'DFT': {
'BASIS_SET_FILE_NAME': 'BASIS_MOLOPT',
'SCF': {
'MAX_SCF': 1,
},
'XC': {
'XC_FUNCTIONAL': {
'_': 'LDA',
},
},
},
'SUBSYS': {
'KIND': {
'_': 'DEFAULT',
'BASIS_SET': 'DZVP-MOLOPT-SR-GTH',
'POTENTIAL': 'GTH-LDA',
},
},
},
})
calc.use_parameters(parameters)

# resources
calc.set_max_wallclock_seconds(3*60) # 3 min
calc.set_resources({"num_machines": 1})

# store and submit
calc.store_all()
calc.submit()
print("submitted calculation: PK=%s" % calc.pk)

wait_for_calc(calc)

# check structure preservation
atoms2 = calc.out.output_structure.get_ase()

# zeros should be preserved exactly
if np.all(atoms2.positions[0] == 0.0):
print("OK, zeros in structure were preserved exactly")
else:
print("ERROR!")
print("Zeros in structure changed: ", atoms2.positions[0])
sys.exit(3)

# other values should be preserved with epsilon precision
dist2 = atoms2.get_distance(0, 1)
if abs(dist2 - dist) < epsilon:
print("OK, structure preserved with %.1e Angstrom precision" % epsilon)
else:
print("ERROR!")
print("Structure changed by %e Angstrom" % abs(dist - dist2))
sys.exit(3)

sys.exit(0)

# EOF

0 comments on commit 96f06d2

Please sign in to comment.