created example for MD reftraj, working

aiidateam · Feb 11, 2024 · 1867f56 · 1867f56
1 parent e763505
commit 1867f56
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 10 deletions.
diff --git a/aiida_cp2k/calculations/__init__.py b/aiida_cp2k/calculations/__init__.py
@@ -6,6 +6,7 @@
 ###############################################################################
 """AiiDA-CP2K input plugin."""
 
+import numpy as np
 from operator import add
 
 from aiida.common import CalcInfo, CodeInfo, InputValidationError
@@ -446,7 +447,7 @@ def _atoms_to_xyz(atoms):
     return xyz
 
 
-def _trajectory_to_xyz_and_cell(tarjectory, atoms):
+def _trajectory_to_xyz_and_cell(trajectory, atoms):
     """Converts postions and cell from a TrajectoryData  to string, taking care of element tags from ASE atoms.
 
     :param atoms: ASE Atoms instance
@@ -457,14 +458,15 @@ def _trajectory_to_xyz_and_cell(tarjectory, atoms):
     xyz = ""
     elem_symbols = kind_names(atoms)
 
-    for step in trajectory.get_array("positions"):
+    for (i,step) in enumerate(trajectory.get_array("positions")):
         elem_coords = [f"{p[0]:25.16f} {p[1]:25.16f} {p[2]:25.16f}" for p in step]
-        xyz += f"{len(elem_coords)}\n\n"
+        xyz += f"{len(elem_coords)}\n"
+        xyz += f"i = {i+1} , time = {(i+1)*0.5} \n"
         xyz += "\n".join(map(add, elem_symbols, elem_coords))
         xyz += "\n"
     if "cells" in trajectory.get_arraynames():
         cell_vecs = [
-            f"{i+1} {(i+1)*0.5:6.3f} {p[0][0]:25.16f} {p[0][1]:25.16f} {p[0][2]:25.16f} {p[1][0]:25.16f} {p[1][1]:25.16f} {p[1][2]:25.16f} {p[2][0]:25.16f} {p[2][1]:25.16f} {p[2][2]:25.16f}"
+            f"{i+1} {(i+1)*0.5:6.3f} {p[0][0]:25.16f} {p[0][1]:25.16f} {p[0][2]:25.16f} {p[1][0]:25.16f} {p[1][1]:25.16f} {p[1][2]:25.16f} {p[2][0]:25.16f} {p[2][1]:25.16f} {p[2][2]:25.16f} {np.dot(p[0],np.cross(p[1],p[2]))}"
             for (i, p) in enumerate(trajectory.get_array("cells"))
         ]
         cell += "\n".join(cell_vecs)

diff --git a/examples/single_calculations/example_dft_md_reftraj.py b/examples/single_calculations/example_dft_md_reftraj.py
@@ -9,8 +9,9 @@
 import os
 import sys
 
-import ase.io
 import click
+import numpy as np
+import ase.io
 from aiida.common import NotExistent
 from aiida.engine import run
 from aiida.orm import Dict, SinglefileData, load_code
@@ -30,12 +31,12 @@ def example_dft_md_reftraj(cp2k_code):
 
     # Structure.
     structure = StructureData(
-        ase=Atoms('HH', positions=[[3, 3, 3], [3.75, 3, 3]]))
+        ase=ase.io.read(os.path.join(thisdir, "..", "files", "h2.xyz"))
     )
 
     # Trajectory.
-    positions = np.array([[[3,3,3],[3.7,3,3]],[[3,3,3],[3.75,3,3]],[[3,3,3],[3.73,3,3]]])
-    cells = np.array([[[6,0,0],[0,6,0],[0,0,6]],[[6.1,0,0],[0,6.2,0],[0,0,6.3]],[[6,0,0],[0,6.1,0],[0,0,6]]])
+    positions = np.array([[[2,2,2.73],[2,2,2.]],[[2,2,2.74],[2,2,2.]],[[2,2,2.75],[2,2,2.]]])
+    cells = np.array([[[4,0,0],[0,4,0],[0,0,4.75]],[[4.4,0,0],[0,4.2,0],[0,0,4.76]],[[4,0,0],[0,4.1,0],[0,0,4.75]]])
     symbols=['H','H']
     trajectory = TrajectoryData()
     trajectory.set_trajectory(symbols, positions, cells=cells)
@@ -56,7 +57,7 @@ def example_dft_md_reftraj(cp2k_code):
             "GLOBAL": {
                 "RUN_TYPE": "MD",
                 "PRINT_LEVEL": "LOW",
-                "WALLTIME": 600
+                "WALLTIME": 600,
                 "PROJECT": "aiida",
             },
             "MOTION": {
@@ -162,7 +163,7 @@ def cli(codelabel):
     except NotExistent:
         print(f"The code '{codelabel}' does not exist.")
         sys.exit(1)
-    example_dft(code)
+    example_dft_md_reftraj(code)
 
 
 if __name__ == "__main__":