Index: ase/calculators/lammps.py
===================================================================
--- ase/calculators/lammps.py	(revision 3)
+++ ase/calculators/lammps.py	(revision 4)
@@ -1,156 +1,568 @@
-<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
-<html xmlns="http://www.w3.org/1999/xhtml">
-  
-  
+#!/usr//bin/env python
 
+# lammps.py (2011/02/22)
+# An ASE calculator for the LAMMPS classical MD code available from
+#       http://lammps.sandia.gov/
+# The environment variable LAMMPS_COMMAND must be defined to point to the LAMMPS binary.
+#
+# Copyright (C) 2009 - 2011 Joerg Meyer, joerg.meyer@ch.tum.de
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this file; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+# or see <http://www.gnu.org/licenses/>.
 
-  <head>
-    <title>
-      lammps.py in trunk/ase/calculators
-     – ase
-    </title>
-    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8" />
-        <link rel="search" href="/projects/ase/search" />
-        <link rel="help" href="/projects/ase/wiki/TracGuide" />
-        <link rel="alternate" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?format=txt" type="text/plain" title="Plain Text" /><link rel="alternate" href="/projects/ase/export/2073/trunk/ase/calculators/lammps.py" type="text/x-python; charset=iso-8859-15" title="Original Format" />
-        <link rel="start" href="/projects/ase/wiki" />
-        <link rel="stylesheet" href="/projects/ase/chrome/common/css/trac.css" type="text/css" /><link rel="stylesheet" href="/projects/ase/chrome/common/css/code.css" type="text/css" /><link rel="stylesheet" href="/projects/ase/chrome/common/css/browser.css" type="text/css" />
-        <link rel="prev" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?rev=2022" title="Revision 2022" />
-        <link rel="shortcut icon" href="/projects/ase/chrome/common/trac.ico" type="image/x-icon" />
-        <link rel="icon" href="/projects/ase/chrome/common/trac.ico" type="image/x-icon" />
-      <link type="application/opensearchdescription+xml" rel="search" href="/projects/ase/search/opensearch" title="Search ase" />
-    <script type="text/javascript" src="/projects/ase/chrome/common/js/jquery.js"></script><script type="text/javascript" src="/projects/ase/chrome/common/js/babel.js"></script><script type="text/javascript" src="/projects/ase/chrome/common/js/trac.js"></script><script type="text/javascript" src="/projects/ase/chrome/common/js/search.js"></script>
-    <!--[if lt IE 7]>
-    <script type="text/javascript" src="/projects/ase/chrome/common/js/ie_pre7_hacks.js"></script>
-    <![endif]-->
-    <script type="text/javascript" src="/projects/ase/chrome/common/js/folding.js"></script><script type="text/javascript">
-      jQuery(document).ready(function($) {
-        $(".trac-toggledeleted").show().click(function() {
-                  $(this).siblings().find(".trac-deleted").toggle();
-                  return false;
-        }).click();
-        $("#jumploc input").hide();
-        $("#jumploc select").change(function () {
-          this.parentNode.parentNode.submit();
-        });
-          $('#preview table.code').enableCollapsibleColumns($('#preview table.code thead th.content'));
-      });
-    </script>
-  </head>
-  <body>
-    <div id="banner">
-      <div id="header">
-        <a id="logo" href="https://trac.fysik.dtu.dk/projects/ase"><img src="https://svn.fysik.dtu.dk/projects/ase/trunk/doc/_static/ase_trac.png" alt="" /></a>
-      </div>
-      <form id="search" action="/projects/ase/search" method="get">
-        <div>
-          <label for="proj-search">Search:</label>
-          <input type="text" id="proj-search" name="q" size="18" value="" />
-          <input type="submit" value="Search" />
-        </div>
-      </form>
-      <div id="metanav" class="nav">
-    <ul>
-      <li class="first"><a href="/projects/ase/login">Login</a></li><li><a href="/projects/ase/wiki/TracGuide">Help/Guide</a></li><li><a href="/projects/ase/about">About Trac</a></li><li class="last"><a href="/projects/ase/prefs">Preferences</a></li>
-    </ul>
-  </div>
-    </div>
-    <div id="mainnav" class="nav">
-    <ul>
-      <li class="first"><a href="/projects/ase/wiki">Wiki</a></li><li><a href="/projects/ase/timeline">Timeline</a></li><li><a href="/projects/ase/roadmap">Roadmap</a></li><li class="active"><a href="/projects/ase/browser">Browse Source</a></li><li><a href="/projects/ase/report">View Tickets</a></li><li class="last"><a href="/projects/ase/search">Search</a></li>
-    </ul>
-  </div>
-    <div id="main">
-      <div id="ctxtnav" class="nav">
-        <h2>Context Navigation</h2>
-          <ul>
-              <li class="first"><span>&larr; <a class="prev" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?rev=2022" title="Revision 2022">Previous Revision</a></span></li><li><span class="missing">Next Revision &rarr;</span></li><li><a href="/projects/ase/browser/trunk/ase/calculators/lammps.py?annotate=blame&amp;rev=2023" title="Annotate each line with the last changed revision (this can be time consuming...)">Annotate</a></li><li class="last"><a href="/projects/ase/log/trunk/ase/calculators/lammps.py">Revision Log</a></li>
-          </ul>
-        <hr />
-      </div>
-    <div id="content" class="browser">
-          <h1>
-<a class="pathentry first" href="/projects/ase/browser?order=name" title="Go to repository root">source:</a>
-<a class="pathentry" href="/projects/ase/browser/trunk?order=name" title="View trunk">trunk</a><span class="pathentry sep">/</span><a class="pathentry" href="/projects/ase/browser/trunk/ase?order=name" title="View ase">ase</a><span class="pathentry sep">/</span><a class="pathentry" href="/projects/ase/browser/trunk/ase/calculators?order=name" title="View calculators">calculators</a><span class="pathentry sep">/</span><a class="pathentry" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?order=name" title="View lammps.py">lammps.py</a>
-<span class="pathentry sep">@</span>
-  <a class="pathentry" href="/projects/ase/changeset/2023" title="View changeset 2023">2023</a>
-<br style="clear: both" />
-</h1>
-        <div id="jumprev">
-          <form action="" method="get">
-            <div>
-              <label for="rev">
-                View revision:</label>
-              <input type="text" id="rev" name="rev" size="6" />
-            </div>
-          </form>
-        </div>
-        <div id="jumploc">
-          <form action="" method="get">
-            <div class="buttons">
-              <label for="preselected">Visit:</label>
-              <select id="preselected" name="preselected">
-                <option selected="selected"></option>
-                <optgroup label="branches">
-                  <option value="/projects/ase/browser/trunk">trunk</option><option value="/projects/ase/browser/branches/dacapo">branches/dacapo</option><option value="/projects/ase/browser/branches/elements">branches/elements</option><option value="/projects/ase/browser/branches/usertags">branches/usertags</option>
-                </optgroup><optgroup label="tags">
-                  <option value="/projects/ase/browser/tags/3.0.0?rev=657">tags/3.0.0</option><option value="/projects/ase/browser/tags/3.1.0?rev=846">tags/3.1.0</option><option value="/projects/ase/browser/tags/3.2.0?rev=1121">tags/3.2.0</option><option value="/projects/ase/browser/tags/3.3.0?rev=1376">tags/3.3.0</option><option value="/projects/ase/browser/tags/3.3.1?rev=1389">tags/3.3.1</option><option value="/projects/ase/browser/tags/3.4.0?rev=1574">tags/3.4.0</option><option value="/projects/ase/browser/tags/3.4.1?rev=1765">tags/3.4.1</option>
-                </optgroup>
-              </select>
-              <input type="submit" value="Go!" title="Jump to the chosen preselected path" />
-            </div>
-          </form>
-        </div>
-      <table id="info" summary="Revision info">
-        <tr>
-          <th scope="col">Revision <a href="/projects/ase/changeset/2023">2023</a>,
-            <span title="23511 bytes">23.0 KB</span>
-            checked in by jensj, <a class="timeline" href="/projects/ase/timeline?from=2011-02-23T08%3A39%3A36%2B01%3A00&amp;precision=second" title="2011-02-23T08:39:36+01:00 in Timeline">4 weeks</a> ago
-            (<a href="/projects/ase/changeset/2023/trunk/ase/calculators/lammps.py">diff</a>)</th>
-        </tr>
-        <tr>
-          <td class="message searchable">
-              <p>
-Python 2.4 compatibility.<br />
-</p>
-          </td>
-        </tr>
-      </table>
-      <div id="preview" class="searchable">
+import os
+import shutil
+import time
+import math
+import decimal as dec
+import numpy as np
+from ase import Atoms
+from ase.parallel import paropen
+
+
+class LAMMPS:
+
+    def __init__(self, label="lammps", dir="LAMMPS", parameters={}, files=[], always_triclinic=False):
+        self.label = label
+        self.dir = dir
+        self.parameters = parameters
+        self.files = files
+        self.always_triclinic = always_triclinic
+        self.calls = 0
+        self.potential_energy = None
+        self.forces = None
+
+        if os.path.isdir(self.dir):
+            if (os.listdir(self.dir) == []):
+                os.rmdir(self.dir)
+            else:
+                os.rename(self.dir, self.dir + ".bak-" + time.strftime("%Y%m%d-%H%M%S"))
+        os.mkdir(self.dir, 0755)
         
-  <table class="code"><thead><tr><th class="lineno" title="Line numbers">Line</th><th class="content"> </th></tr></thead><tbody><tr><th id="L1"><a href="#L1">1</a></th><td>#!/usr//bin/env python</td></tr><tr><th id="L2"><a href="#L2">2</a></th><td></td></tr><tr><th id="L3"><a href="#L3">3</a></th><td># lammps.py (2011/02/22)</td></tr><tr><th id="L4"><a href="#L4">4</a></th><td># An ASE calculator for the LAMMPS classical MD code available from</td></tr><tr><th id="L5"><a href="#L5">5</a></th><td>#       http://lammps.sandia.gov/</td></tr><tr><th id="L6"><a href="#L6">6</a></th><td># The environment variable LAMMPS_COMMAND must be defined to point to the LAMMPS binary.</td></tr><tr><th id="L7"><a href="#L7">7</a></th><td>#</td></tr><tr><th id="L8"><a href="#L8">8</a></th><td># Copyright (C) 2009 - 2011 Joerg Meyer, joerg.meyer@ch.tum.de</td></tr><tr><th id="L9"><a href="#L9">9</a></th><td>#</td></tr><tr><th id="L10"><a href="#L10">10</a></th><td># This library is free software; you can redistribute it and/or</td></tr><tr><th id="L11"><a href="#L11">11</a></th><td># modify it under the terms of the GNU Lesser General Public</td></tr><tr><th id="L12"><a href="#L12">12</a></th><td># License as published by the Free Software Foundation; either</td></tr><tr><th id="L13"><a href="#L13">13</a></th><td># version 2.1 of the License, or (at your option) any later version.</td></tr><tr><th id="L14"><a href="#L14">14</a></th><td>#</td></tr><tr><th id="L15"><a href="#L15">15</a></th><td># This library is distributed in the hope that it will be useful,</td></tr><tr><th id="L16"><a href="#L16">16</a></th><td># but WITHOUT ANY WARRANTY; without even the implied warranty of</td></tr><tr><th id="L17"><a href="#L17">17</a></th><td># MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU</td></tr><tr><th id="L18"><a href="#L18">18</a></th><td># Lesser General Public License for more details.</td></tr><tr><th id="L19"><a href="#L19">19</a></th><td>#</td></tr><tr><th id="L20"><a href="#L20">20</a></th><td># You should have received a copy of the GNU Lesser General Public</td></tr><tr><th id="L21"><a href="#L21">21</a></th><td># License along with this file; if not, write to the Free Software</td></tr><tr><th id="L22"><a href="#L22">22</a></th><td># Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA</td></tr><tr><th id="L23"><a href="#L23">23</a></th><td># or see &lt;http://www.gnu.org/licenses/&gt;.</td></tr><tr><th id="L24"><a href="#L24">24</a></th><td></td></tr><tr><th id="L25"><a href="#L25">25</a></th><td>import os</td></tr><tr><th id="L26"><a href="#L26">26</a></th><td>import shutil</td></tr><tr><th id="L27"><a href="#L27">27</a></th><td>import time</td></tr><tr><th id="L28"><a href="#L28">28</a></th><td>import math</td></tr><tr><th id="L29"><a href="#L29">29</a></th><td>import decimal as dec</td></tr><tr><th id="L30"><a href="#L30">30</a></th><td>import numpy as np</td></tr><tr><th id="L31"><a href="#L31">31</a></th><td>from ase import Atoms</td></tr><tr><th id="L32"><a href="#L32">32</a></th><td>from ase.parallel import paropen</td></tr><tr><th id="L33"><a href="#L33">33</a></th><td></td></tr><tr><th id="L34"><a href="#L34">34</a></th><td></td></tr><tr><th id="L35"><a href="#L35">35</a></th><td>class LAMMPS:</td></tr><tr><th id="L36"><a href="#L36">36</a></th><td></td></tr><tr><th id="L37"><a href="#L37">37</a></th><td>    def __init__(self, label="lammps", dir="LAMMPS", parameters={}, files=[], always_triclinic=False):</td></tr><tr><th id="L38"><a href="#L38">38</a></th><td>        self.label = label</td></tr><tr><th id="L39"><a href="#L39">39</a></th><td>        self.dir = dir</td></tr><tr><th id="L40"><a href="#L40">40</a></th><td>        self.parameters = parameters</td></tr><tr><th id="L41"><a href="#L41">41</a></th><td>        self.files = files</td></tr><tr><th id="L42"><a href="#L42">42</a></th><td>        self.always_triclinic = always_triclinic</td></tr><tr><th id="L43"><a href="#L43">43</a></th><td>        self.calls = 0</td></tr><tr><th id="L44"><a href="#L44">44</a></th><td>        self.potential_energy = None</td></tr><tr><th id="L45"><a href="#L45">45</a></th><td>        self.forces = None</td></tr><tr><th id="L46"><a href="#L46">46</a></th><td></td></tr><tr><th id="L47"><a href="#L47">47</a></th><td>        if os.path.isdir(self.dir):</td></tr><tr><th id="L48"><a href="#L48">48</a></th><td>            if (os.listdir(self.dir) == []):</td></tr><tr><th id="L49"><a href="#L49">49</a></th><td>                os.rmdir(self.dir)</td></tr><tr><th id="L50"><a href="#L50">50</a></th><td>            else:</td></tr><tr><th id="L51"><a href="#L51">51</a></th><td>                os.rename(self.dir, self.dir + ".bak-" + time.strftime("%Y%m%d-%H%M%S"))</td></tr><tr><th id="L52"><a href="#L52">52</a></th><td>        os.mkdir(self.dir, 0755)</td></tr><tr><th id="L53"><a href="#L53">53</a></th><td>        </td></tr><tr><th id="L54"><a href="#L54">54</a></th><td>        for f in files:</td></tr><tr><th id="L55"><a href="#L55">55</a></th><td>            shutil.copy(f, os.path.join(self.dir, f))</td></tr><tr><th id="L56"><a href="#L56">56</a></th><td></td></tr><tr><th id="L57"><a href="#L57">57</a></th><td>    def clean(self, force=False):</td></tr><tr><th id="L58"><a href="#L58">58</a></th><td>        if os.path.isdir(self.dir):</td></tr><tr><th id="L59"><a href="#L59">59</a></th><td>            files_in_dir = os.listdir(self.dir)</td></tr><tr><th id="L60"><a href="#L60">60</a></th><td>            files_copied = self.files</td></tr><tr><th id="L61"><a href="#L61">61</a></th><td>            if (len(files_in_dir) == 0):</td></tr><tr><th id="L62"><a href="#L62">62</a></th><td>                os.rmdir(self.dir)</td></tr><tr><th id="L63"><a href="#L63">63</a></th><td>            elif (set(files_in_dir).issubset(set(files_copied)) or force):</td></tr><tr><th id="L64"><a href="#L64">64</a></th><td>                shutil.rmtree(self.dir)</td></tr><tr><th id="L65"><a href="#L65">65</a></th><td></td></tr><tr><th id="L66"><a href="#L66">66</a></th><td>    def get_potential_energy(self, atoms):</td></tr><tr><th id="L67"><a href="#L67">67</a></th><td>        self.update(atoms)</td></tr><tr><th id="L68"><a href="#L68">68</a></th><td>        return self.potential_energy</td></tr><tr><th id="L69"><a href="#L69">69</a></th><td></td></tr><tr><th id="L70"><a href="#L70">70</a></th><td>    def get_forces(self, atoms):</td></tr><tr><th id="L71"><a href="#L71">71</a></th><td>        self.update(atoms)</td></tr><tr><th id="L72"><a href="#L72">72</a></th><td>        return self.forces.copy()</td></tr><tr><th id="L73"><a href="#L73">73</a></th><td></td></tr><tr><th id="L74"><a href="#L74">74</a></th><td>    def get_stress(self, atoms):</td></tr><tr><th id="L75"><a href="#L75">75</a></th><td>        raise NotImplementedError</td></tr><tr><th id="L76"><a href="#L76">76</a></th><td></td></tr><tr><th id="L77"><a href="#L77">77</a></th><td>    def update(self, atoms):</td></tr><tr><th id="L78"><a href="#L78">78</a></th><td>        # TODO: check if (re-)calculation is necessary</td></tr><tr><th id="L79"><a href="#L79">79</a></th><td>        self.calculate(atoms)</td></tr><tr><th id="L80"><a href="#L80">80</a></th><td></td></tr><tr><th id="L81"><a href="#L81">81</a></th><td>    def calculate(self, atoms):</td></tr><tr><th id="L82"><a href="#L82">82</a></th><td>        self.atoms = atoms.copy()</td></tr><tr><th id="L83"><a href="#L83">83</a></th><td>        self.run()</td></tr><tr><th id="L84"><a href="#L84">84</a></th><td></td></tr><tr><th id="L85"><a href="#L85">85</a></th><td>    def run(self):</td></tr><tr><th id="L86"><a href="#L86">86</a></th><td>        """Method which explicitely runs LAMMPS."""</td></tr><tr><th id="L87"><a href="#L87">87</a></th><td></td></tr><tr><th id="L88"><a href="#L88">88</a></th><td>        self.calls += 1</td></tr><tr><th id="L89"><a href="#L89">89</a></th><td></td></tr><tr><th id="L90"><a href="#L90">90</a></th><td>        # set LAMMPS command from environment variable</td></tr><tr><th id="L91"><a href="#L91">91</a></th><td>        if os.environ.has_key('LAMMPS_COMMAND'):</td></tr><tr><th id="L92"><a href="#L92">92</a></th><td>            lammps_command = os.path.abspath(os.environ['LAMMPS_COMMAND'])</td></tr><tr><th id="L93"><a href="#L93">93</a></th><td>        else:</td></tr><tr><th id="L94"><a href="#L94">94</a></th><td>            self.clean()</td></tr><tr><th id="L95"><a href="#L95">95</a></th><td>            raise RuntimeError('Please set LAMMPS_COMMAND environment variable')</td></tr><tr><th id="L96"><a href="#L96">96</a></th><td>        if os.environ.has_key('LAMMPS_OPTIONS'):</td></tr><tr><th id="L97"><a href="#L97">97</a></th><td>            lammps_options = os.environ['LAMMPS_OPTIONS']</td></tr><tr><th id="L98"><a href="#L98">98</a></th><td>        else:</td></tr><tr><th id="L99"><a href="#L99">99</a></th><td>            lammps_options = "-echo log -screen none"</td></tr><tr><th id="L100"><a href="#L100">100</a></th><td></td></tr><tr><th id="L101"><a href="#L101">101</a></th><td>        # change into subdirectory for LAMMPS calculations</td></tr><tr><th id="L102"><a href="#L102">102</a></th><td>        cwd = os.getcwd()</td></tr><tr><th id="L103"><a href="#L103">103</a></th><td>        os.chdir(self.dir)</td></tr><tr><th id="L104"><a href="#L104">104</a></th><td></td></tr><tr><th id="L105"><a href="#L105">105</a></th><td>        # setup file names for LAMMPS calculation</td></tr><tr><th id="L106"><a href="#L106">106</a></th><td>        label = "%s-%06d" % (self.label, self.calls)</td></tr><tr><th id="L107"><a href="#L107">107</a></th><td>        lammps_data = "data." + label</td></tr><tr><th id="L108"><a href="#L108">108</a></th><td>        lammps_in = "in." + label</td></tr><tr><th id="L109"><a href="#L109">109</a></th><td>        lammps_trj = label + ".lammpstrj"</td></tr><tr><th id="L110"><a href="#L110">110</a></th><td>        lammps_log = label + ".log"</td></tr><tr><th id="L111"><a href="#L111">111</a></th><td></td></tr><tr><th id="L112"><a href="#L112">112</a></th><td>        # write LAMMPS input files</td></tr><tr><th id="L113"><a href="#L113">113</a></th><td>        self.write_lammps_data(lammps_data=lammps_data)</td></tr><tr><th id="L114"><a href="#L114">114</a></th><td>        self.write_lammps_in(lammps_in=lammps_in, lammps_data=lammps_data, lammps_trj=lammps_trj, parameters=self.parameters)</td></tr><tr><th id="L115"><a href="#L115">115</a></th><td></td></tr><tr><th id="L116"><a href="#L116">116</a></th><td>        # run LAMMPS</td></tr><tr><th id="L117"><a href="#L117">117</a></th><td>        # TODO: check for successful completion (based on output files?!) of LAMMPS call...</td></tr><tr><th id="L118"><a href="#L118">118</a></th><td>        exitcode = os.system("%s %s -in %s -log %s" % (lammps_command, lammps_options, lammps_in, lammps_log))</td></tr><tr><th id="L119"><a href="#L119">119</a></th><td>        if exitcode != 0:</td></tr><tr><th id="L120"><a href="#L120">120</a></th><td>            cwd = os.getcwd()</td></tr><tr><th id="L121"><a href="#L121">121</a></th><td>            raise RuntimeError("LAMMPS exited in %s with exit code: %d.  " % (cwd,exitcode))</td></tr><tr><th id="L122"><a href="#L122">122</a></th><td></td></tr><tr><th id="L123"><a href="#L123">123</a></th><td>        # read LAMMPS output files</td></tr><tr><th id="L124"><a href="#L124">124</a></th><td>        self.read_lammps_log(lammps_log=lammps_log)</td></tr><tr><th id="L125"><a href="#L125">125</a></th><td>        self.read_lammps_trj(lammps_trj=lammps_trj)</td></tr><tr><th id="L126"><a href="#L126">126</a></th><td></td></tr><tr><th id="L127"><a href="#L127">127</a></th><td>        # change back to previous working directory</td></tr><tr><th id="L128"><a href="#L128">128</a></th><td>        os.chdir(cwd)</td></tr><tr><th id="L129"><a href="#L129">129</a></th><td></td></tr><tr><th id="L130"><a href="#L130">130</a></th><td>    def write_lammps_data(self, lammps_data=None):</td></tr><tr><th id="L131"><a href="#L131">131</a></th><td>        """Method which writes a LAMMPS data file with atomic structure."""</td></tr><tr><th id="L132"><a href="#L132">132</a></th><td>        if (lammps_data == None):</td></tr><tr><th id="L133"><a href="#L133">133</a></th><td>            lammps_data = "data." + self.label</td></tr><tr><th id="L134"><a href="#L134">134</a></th><td>        write_lammps(lammps_data, self.atoms, force_skew=self.always_triclinic)</td></tr><tr><th id="L135"><a href="#L135">135</a></th><td></td></tr><tr><th id="L136"><a href="#L136">136</a></th><td>    def write_lammps_in(self, lammps_in=None, lammps_data=None, lammps_trj=None, parameters={}):</td></tr><tr><th id="L137"><a href="#L137">137</a></th><td>        """Method which writes a LAMMPS in file with run parameters and settings."""</td></tr><tr><th id="L138"><a href="#L138">138</a></th><td>        if (lammps_in == None):</td></tr><tr><th id="L139"><a href="#L139">139</a></th><td>            lammps_in = "in." + self.label</td></tr><tr><th id="L140"><a href="#L140">140</a></th><td>        if (lammps_data == None):</td></tr><tr><th id="L141"><a href="#L141">141</a></th><td>            lammps_data = "data." + self.label</td></tr><tr><th id="L142"><a href="#L142">142</a></th><td>        if (lammps_trj == None):</td></tr><tr><th id="L143"><a href="#L143">143</a></th><td>            lammps_trj = self.label + ".lammpstrj"</td></tr><tr><th id="L144"><a href="#L144">144</a></th><td></td></tr><tr><th id="L145"><a href="#L145">145</a></th><td>        f = paropen(lammps_in, 'w')</td></tr><tr><th id="L146"><a href="#L146">146</a></th><td>        f.write("# " + f.name + " (written by ASE) \n")</td></tr><tr><th id="L147"><a href="#L147">147</a></th><td>        f.write("\n")</td></tr><tr><th id="L148"><a href="#L148">148</a></th><td></td></tr><tr><th id="L149"><a href="#L149">149</a></th><td>        f.write("### variables \n")</td></tr><tr><th id="L150"><a href="#L150">150</a></th><td>        f.write("variable data_file index \"%s\" \n" % lammps_data)</td></tr><tr><th id="L151"><a href="#L151">151</a></th><td>        f.write("variable dump_file index \"%s\" \n" % lammps_trj)</td></tr><tr><th id="L152"><a href="#L152">152</a></th><td>        f.write("\n\n")</td></tr><tr><th id="L153"><a href="#L153">153</a></th><td></td></tr><tr><th id="L154"><a href="#L154">154</a></th><td>        pbc = self.atoms.get_pbc()</td></tr><tr><th id="L155"><a href="#L155">155</a></th><td>        f.write("### simulation box \n")</td></tr><tr><th id="L156"><a href="#L156">156</a></th><td>        f.write("units metal \n")</td></tr><tr><th id="L157"><a href="#L157">157</a></th><td>        if ("boundary" in parameters):</td></tr><tr><th id="L158"><a href="#L158">158</a></th><td>            f.write("boundary %s \n" % parameters["boundary"])</td></tr><tr><th id="L159"><a href="#L159">159</a></th><td>        else:</td></tr><tr><th id="L160"><a href="#L160">160</a></th><td>            f.write("boundary %c %c %c \n" % tuple('sp'[x] for x in pbc))</td></tr><tr><th id="L161"><a href="#L161">161</a></th><td>        f.write("atom_modify sort 0 0.0 \n")</td></tr><tr><th id="L162"><a href="#L162">162</a></th><td>        for key in ("neighbor" ,"newton"):</td></tr><tr><th id="L163"><a href="#L163">163</a></th><td>            if key in parameters:</td></tr><tr><th id="L164"><a href="#L164">164</a></th><td>                f.write("%s %s \n" % (key, parameters[key]))</td></tr><tr><th id="L165"><a href="#L165">165</a></th><td>        f.write("\n")</td></tr><tr><th id="L166"><a href="#L166">166</a></th><td>        f.write("read_data ${data_file} \n")</td></tr><tr><th id="L167"><a href="#L167">167</a></th><td>        f.write("\n\n")</td></tr><tr><th id="L168"><a href="#L168">168</a></th><td></td></tr><tr><th id="L169"><a href="#L169">169</a></th><td>        f.write("### interactions \n")</td></tr><tr><th id="L170"><a href="#L170">170</a></th><td>        if ( ("pair_style" in parameters) and ("pair_coeff" in parameters) ):</td></tr><tr><th id="L171"><a href="#L171">171</a></th><td>            pair_style = parameters["pair_style"]</td></tr><tr><th id="L172"><a href="#L172">172</a></th><td>            f.write("pair_style %s \n" % pair_style)</td></tr><tr><th id="L173"><a href="#L173">173</a></th><td>            for pair_coeff in parameters["pair_coeff"]:</td></tr><tr><th id="L174"><a href="#L174">174</a></th><td>                f.write("pair_coeff %s \n" % pair_coeff)</td></tr><tr><th id="L175"><a href="#L175">175</a></th><td>            if "mass" in parameters:</td></tr><tr><th id="L176"><a href="#L176">176</a></th><td>                for mass in parameters["mass"]:</td></tr><tr><th id="L177"><a href="#L177">177</a></th><td>                    f.write("mass %s \n" % mass)</td></tr><tr><th id="L178"><a href="#L178">178</a></th><td>        else:</td></tr><tr><th id="L179"><a href="#L179">179</a></th><td>            # simple default parameters that should always make the LAMMPS calculation run</td></tr><tr><th id="L180"><a href="#L180">180</a></th><td>            #    pair_style     lj/cut 2.5</td></tr><tr><th id="L181"><a href="#L181">181</a></th><td>            #    pair_coeff     * * 1 1</td></tr><tr><th id="L182"><a href="#L182">182</a></th><td>            #    mass           * 1.0        else:</td></tr><tr><th id="L183"><a href="#L183">183</a></th><td>            f.write("pair_style lj/cut 2.5 \n")</td></tr><tr><th id="L184"><a href="#L184">184</a></th><td>            f.write("pair_coeff * * 1 1 \n")</td></tr><tr><th id="L185"><a href="#L185">185</a></th><td>            f.write("mass * 1.0 \n")</td></tr><tr><th id="L186"><a href="#L186">186</a></th><td>        f.write("\n")</td></tr><tr><th id="L187"><a href="#L187">187</a></th><td></td></tr><tr><th id="L188"><a href="#L188">188</a></th><td>        f.write("### run \n")</td></tr><tr><th id="L189"><a href="#L189">189</a></th><td>        f.write("fix fix_nve all nve \n")</td></tr><tr><th id="L190"><a href="#L190">190</a></th><td>        f.write("\n")</td></tr><tr><th id="L191"><a href="#L191">191</a></th><td>        f.write("dump dump_all all custom 1 ${dump_file} id type x y z vx vy vz fx fy fz \n")</td></tr><tr><th id="L192"><a href="#L192">192</a></th><td>        f.write("\n")</td></tr><tr><th id="L193"><a href="#L193">193</a></th><td>        f.write("thermo_style custom step temp ke pe etotal cpu \n")</td></tr><tr><th id="L194"><a href="#L194">194</a></th><td>        f.write("thermo_modify format 1 %4d format 2 %9.3f format 3 %20.6f format 4 %20.6f format 5 %20.6f format 6 %9.3f \n")</td></tr><tr><th id="L195"><a href="#L195">195</a></th><td>        f.write("thermo 1 \n")</td></tr><tr><th id="L196"><a href="#L196">196</a></th><td>        f.write("\n")</td></tr><tr><th id="L197"><a href="#L197">197</a></th><td>        if ("minimize" in parameters):</td></tr><tr><th id="L198"><a href="#L198">198</a></th><td>            f.write("minimize %s \n" % parameters["minimize"])</td></tr><tr><th id="L199"><a href="#L199">199</a></th><td>        if ("run" in parameters):</td></tr><tr><th id="L200"><a href="#L200">200</a></th><td>            f.write("run %s \n" % parameters["run"])</td></tr><tr><th id="L201"><a href="#L201">201</a></th><td>        if not ( ("minimize" in parameters) or ("run" in parameters) ):</td></tr><tr><th id="L202"><a href="#L202">202</a></th><td>            f.write("run 0 \n")</td></tr><tr><th id="L203"><a href="#L203">203</a></th><td></td></tr><tr><th id="L204"><a href="#L204">204</a></th><td>        f.close()</td></tr><tr><th id="L205"><a href="#L205">205</a></th><td></td></tr><tr><th id="L206"><a href="#L206">206</a></th><td>    def read_lammps_log(self, lammps_log=None):</td></tr><tr><th id="L207"><a href="#L207">207</a></th><td>        """Method which reads a LAMMPS output log file."""</td></tr><tr><th id="L208"><a href="#L208">208</a></th><td>        if (lammps_log == None):</td></tr><tr><th id="L209"><a href="#L209">209</a></th><td>            lammps_log = self.label + ".log"</td></tr><tr><th id="L210"><a href="#L210">210</a></th><td></td></tr><tr><th id="L211"><a href="#L211">211</a></th><td>        epot = 0.0</td></tr><tr><th id="L212"><a href="#L212">212</a></th><td></td></tr><tr><th id="L213"><a href="#L213">213</a></th><td>        f = paropen(lammps_log, 'r')</td></tr><tr><th id="L214"><a href="#L214">214</a></th><td>        while True:</td></tr><tr><th id="L215"><a href="#L215">215</a></th><td>            line = f.readline()</td></tr><tr><th id="L216"><a href="#L216">216</a></th><td>            if not line:</td></tr><tr><th id="L217"><a href="#L217">217</a></th><td>                break</td></tr><tr><th id="L218"><a href="#L218">218</a></th><td>            # get potential energy of first step (if more are done)</td></tr><tr><th id="L219"><a href="#L219">219</a></th><td>            if "PotEng" in line:</td></tr><tr><th id="L220"><a href="#L220">220</a></th><td>                i = line.split().index("PotEng")</td></tr><tr><th id="L221"><a href="#L221">221</a></th><td>                line = f.readline()</td></tr><tr><th id="L222"><a href="#L222">222</a></th><td>                epot = float(line.split()[i])</td></tr><tr><th id="L223"><a href="#L223">223</a></th><td>        f.close()</td></tr><tr><th id="L224"><a href="#L224">224</a></th><td></td></tr><tr><th id="L225"><a href="#L225">225</a></th><td>#        print epot</td></tr><tr><th id="L226"><a href="#L226">226</a></th><td></td></tr><tr><th id="L227"><a href="#L227">227</a></th><td>        self.potential_energy = epot</td></tr><tr><th id="L228"><a href="#L228">228</a></th><td></td></tr><tr><th id="L229"><a href="#L229">229</a></th><td>    def read_lammps_trj(self, lammps_trj=None, set_atoms=False):</td></tr><tr><th id="L230"><a href="#L230">230</a></th><td>        """Method which reads a LAMMPS dump file."""</td></tr><tr><th id="L231"><a href="#L231">231</a></th><td>        if (lammps_trj == None):</td></tr><tr><th id="L232"><a href="#L232">232</a></th><td>            lammps_trj = self.label + ".lammpstrj"</td></tr><tr><th id="L233"><a href="#L233">233</a></th><td></td></tr><tr><th id="L234"><a href="#L234">234</a></th><td>        f = paropen(lammps_trj, 'r')</td></tr><tr><th id="L235"><a href="#L235">235</a></th><td>        while True:</td></tr><tr><th id="L236"><a href="#L236">236</a></th><td>            line = f.readline()</td></tr><tr><th id="L237"><a href="#L237">237</a></th><td></td></tr><tr><th id="L238"><a href="#L238">238</a></th><td>            if not line:</td></tr><tr><th id="L239"><a href="#L239">239</a></th><td>                break</td></tr><tr><th id="L240"><a href="#L240">240</a></th><td></td></tr><tr><th id="L241"><a href="#L241">241</a></th><td>            #TODO: extend to proper dealing with multiple steps in one trajectory file</td></tr><tr><th id="L242"><a href="#L242">242</a></th><td>            if "ITEM: TIMESTEP" in line:</td></tr><tr><th id="L243"><a href="#L243">243</a></th><td>                n_atoms = 0</td></tr><tr><th id="L244"><a href="#L244">244</a></th><td>                lo = [] ; hi = [] ; tilt = []</td></tr><tr><th id="L245"><a href="#L245">245</a></th><td>                id = [] ; type = []</td></tr><tr><th id="L246"><a href="#L246">246</a></th><td>                positions = [] ; velocities = [] ; forces = []</td></tr><tr><th id="L247"><a href="#L247">247</a></th><td></td></tr><tr><th id="L248"><a href="#L248">248</a></th><td>            if "ITEM: NUMBER OF ATOMS" in line:</td></tr><tr><th id="L249"><a href="#L249">249</a></th><td>                line = f.readline()</td></tr><tr><th id="L250"><a href="#L250">250</a></th><td>                n_atoms = int(line.split()[0])</td></tr><tr><th id="L251"><a href="#L251">251</a></th><td>            </td></tr><tr><th id="L252"><a href="#L252">252</a></th><td>            if "ITEM: BOX BOUNDS" in line:</td></tr><tr><th id="L253"><a href="#L253">253</a></th><td>                # save labels behind "ITEM: BOX BOUNDS" in triclinic case (&gt;=lammps-7Jul09)</td></tr><tr><th id="L254"><a href="#L254">254</a></th><td>                tilt_items = line.split()[3:]</td></tr><tr><th id="L255"><a href="#L255">255</a></th><td>                for i in range(3):</td></tr><tr><th id="L256"><a href="#L256">256</a></th><td>                    line = f.readline()</td></tr><tr><th id="L257"><a href="#L257">257</a></th><td>                    fields = line.split()</td></tr><tr><th id="L258"><a href="#L258">258</a></th><td>                    lo.append(float(fields[0]))</td></tr><tr><th id="L259"><a href="#L259">259</a></th><td>                    hi.append(float(fields[1]))</td></tr><tr><th id="L260"><a href="#L260">260</a></th><td>                    if (len(fields) &gt;= 3):</td></tr><tr><th id="L261"><a href="#L261">261</a></th><td>                        tilt.append(float(fields[2]))</td></tr><tr><th id="L262"><a href="#L262">262</a></th><td>            </td></tr><tr><th id="L263"><a href="#L263">263</a></th><td>            if "ITEM: ATOMS" in line:</td></tr><tr><th id="L264"><a href="#L264">264</a></th><td>                # (reliably) identify values by labels behind "ITEM: ATOMS" - requires &gt;=lammps-7Jul09</td></tr><tr><th id="L265"><a href="#L265">265</a></th><td>                # create corresponding index dictionary before iterating over atoms to (hopefully) speed up lookups...</td></tr><tr><th id="L266"><a href="#L266">266</a></th><td>                atom_attributes = {}</td></tr><tr><th id="L267"><a href="#L267">267</a></th><td>                for (i, x) in enumerate(line.split()[2:]):</td></tr><tr><th id="L268"><a href="#L268">268</a></th><td>                    atom_attributes[x] = i</td></tr><tr><th id="L269"><a href="#L269">269</a></th><td>                for n in range(n_atoms):</td></tr><tr><th id="L270"><a href="#L270">270</a></th><td>                    line = f.readline()</td></tr><tr><th id="L271"><a href="#L271">271</a></th><td>                    fields = line.split()</td></tr><tr><th id="L272"><a href="#L272">272</a></th><td>                    id.append( int(fields[atom_attributes['id']]) )</td></tr><tr><th id="L273"><a href="#L273">273</a></th><td>                    type.append( int(fields[atom_attributes['type']]) )</td></tr><tr><th id="L274"><a href="#L274">274</a></th><td>                    positions.append( [ float(fields[atom_attributes[x]]) for x in ['x', 'y', 'z'] ] )</td></tr><tr><th id="L275"><a href="#L275">275</a></th><td>                    velocities.append( [ float(fields[atom_attributes[x]]) for x in ['vx', 'vy', 'vz'] ] )</td></tr><tr><th id="L276"><a href="#L276">276</a></th><td>                    forces.append( [ float(fields[atom_attributes[x]]) for x in ['fx', 'fy', 'fz'] ] )</td></tr><tr><th id="L277"><a href="#L277">277</a></th><td>        f.close()</td></tr><tr><th id="L278"><a href="#L278">278</a></th><td></td></tr><tr><th id="L279"><a href="#L279">279</a></th><td>        # determine cell tilt (triclinic case!)</td></tr><tr><th id="L280"><a href="#L280">280</a></th><td>        if (len(tilt) &gt;= 3):</td></tr><tr><th id="L281"><a href="#L281">281</a></th><td>            # for &gt;=lammps-7Jul09 use labels behind "ITEM: BOX BOUNDS" to assign tilt (vector) elements ...</td></tr><tr><th id="L282"><a href="#L282">282</a></th><td>            if (len(tilt_items) &gt;= 3):</td></tr><tr><th id="L283"><a href="#L283">283</a></th><td>                xy = tilt[tilt_items.index('xy')]</td></tr><tr><th id="L284"><a href="#L284">284</a></th><td>                xz = tilt[tilt_items.index('xz')]</td></tr><tr><th id="L285"><a href="#L285">285</a></th><td>                yz = tilt[tilt_items.index('yz')]</td></tr><tr><th id="L286"><a href="#L286">286</a></th><td>            # ... otherwise assume default order in 3rd column (if the latter was present)</td></tr><tr><th id="L287"><a href="#L287">287</a></th><td>            else:</td></tr><tr><th id="L288"><a href="#L288">288</a></th><td>                xy = tilt[0]</td></tr><tr><th id="L289"><a href="#L289">289</a></th><td>                xz = tilt[1]</td></tr><tr><th id="L290"><a href="#L290">290</a></th><td>                yz = tilt[2]</td></tr><tr><th id="L291"><a href="#L291">291</a></th><td>        else:</td></tr><tr><th id="L292"><a href="#L292">292</a></th><td>            xy = xz = yz = 0</td></tr><tr><th id="L293"><a href="#L293">293</a></th><td>        xhilo = (hi[0] - lo[0]) - xy - xz</td></tr><tr><th id="L294"><a href="#L294">294</a></th><td>        yhilo = (hi[1] - lo[1]) - yz</td></tr><tr><th id="L295"><a href="#L295">295</a></th><td>        zhilo = (hi[2] - lo[2])</td></tr><tr><th id="L296"><a href="#L296">296</a></th><td>        </td></tr><tr><th id="L297"><a href="#L297">297</a></th><td>        # The simulation box bounds are included in each snapshot and if the box is triclinic (non-orthogonal), </td></tr><tr><th id="L298"><a href="#L298">298</a></th><td>        # then the tilt factors are also printed; see the region prism command for a description of tilt factors. </td></tr><tr><th id="L299"><a href="#L299">299</a></th><td>        # For triclinic boxes the box bounds themselves (first 2 quantities on each line) are a true "bounding box" </td></tr><tr><th id="L300"><a href="#L300">300</a></th><td>        # around the simulation domain, which means they include the effect of any tilt.</td></tr><tr><th id="L301"><a href="#L301">301</a></th><td>        # [ http://lammps.sandia.gov/doc/dump.html , lammps-7Jul09 ]</td></tr><tr><th id="L302"><a href="#L302">302</a></th><td>        #</td></tr><tr><th id="L303"><a href="#L303">303</a></th><td>        # This *should* extract the lattice vectors that LAMMPS uses from the true "bounding box" printed in the dump file</td></tr><tr><th id="L304"><a href="#L304">304</a></th><td>        # It might fail in some cases (negative tilts?!) due to the MIN / MAX construction of these box corners:</td></tr><tr><th id="L305"><a href="#L305">305</a></th><td>        #</td></tr><tr><th id="L306"><a href="#L306">306</a></th><td>        #       void Domain::set_global_box() </td></tr><tr><th id="L307"><a href="#L307">307</a></th><td>        #       [...]</td></tr><tr><th id="L308"><a href="#L308">308</a></th><td>        #         if (triclinic) {</td></tr><tr><th id="L309"><a href="#L309">309</a></th><td>        #           [...]</td></tr><tr><th id="L310"><a href="#L310">310</a></th><td>        #           boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy);</td></tr><tr><th id="L311"><a href="#L311">311</a></th><td>        #           boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz);</td></tr><tr><th id="L312"><a href="#L312">312</a></th><td>        #           boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz);</td></tr><tr><th id="L313"><a href="#L313">313</a></th><td>        #           boxlo_bound[2] = boxlo[2];</td></tr><tr><th id="L314"><a href="#L314">314</a></th><td>        #</td></tr><tr><th id="L315"><a href="#L315">315</a></th><td>        #           boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy);</td></tr><tr><th id="L316"><a href="#L316">316</a></th><td>        #           boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz);</td></tr><tr><th id="L317"><a href="#L317">317</a></th><td>        #           boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);</td></tr><tr><th id="L318"><a href="#L318">318</a></th><td>        #           boxhi_bound[2] = boxhi[2];</td></tr><tr><th id="L319"><a href="#L319">319</a></th><td>        #         }</td></tr><tr><th id="L320"><a href="#L320">320</a></th><td>        # [ lammps-7Jul09/src/domain.cpp ]</td></tr><tr><th id="L321"><a href="#L321">321</a></th><td>        #</td></tr><tr><th id="L322"><a href="#L322">322</a></th><td>        cell = [[xhilo,0,0],[xy,yhilo,0],[xz,yz,zhilo]]</td></tr><tr><th id="L323"><a href="#L323">323</a></th><td></td></tr><tr><th id="L324"><a href="#L324">324</a></th><td>#       print n_atoms</td></tr><tr><th id="L325"><a href="#L325">325</a></th><td>#       print lo</td></tr><tr><th id="L326"><a href="#L326">326</a></th><td>#       print hi</td></tr><tr><th id="L327"><a href="#L327">327</a></th><td>#       print id</td></tr><tr><th id="L328"><a href="#L328">328</a></th><td>#       print type</td></tr><tr><th id="L329"><a href="#L329">329</a></th><td>#       print positions</td></tr><tr><th id="L330"><a href="#L330">330</a></th><td>#       print velocities</td></tr><tr><th id="L331"><a href="#L331">331</a></th><td>#       print forces</td></tr><tr><th id="L332"><a href="#L332">332</a></th><td></td></tr><tr><th id="L333"><a href="#L333">333</a></th><td>        # assume that LAMMPS does not reorder atoms internally</td></tr><tr><th id="L334"><a href="#L334">334</a></th><td>#        print np.array(positions)</td></tr><tr><th id="L335"><a href="#L335">335</a></th><td>#        print self.atoms.get_positions()</td></tr><tr><th id="L336"><a href="#L336">336</a></th><td>#        print np.array(positions) - self.atoms.get_positions()</td></tr><tr><th id="L337"><a href="#L337">337</a></th><td></td></tr><tr><th id="L338"><a href="#L338">338</a></th><td>        cell_atoms = np.array(cell)</td></tr><tr><th id="L339"><a href="#L339">339</a></th><td>        type_atoms = np.array(type)</td></tr><tr><th id="L340"><a href="#L340">340</a></th><td>#      velocities_atoms = np.array(velocities)</td></tr><tr><th id="L341"><a href="#L341">341</a></th><td>        positions_atoms = np.array(positions)</td></tr><tr><th id="L342"><a href="#L342">342</a></th><td>        forces_atoms = np.array(forces)</td></tr><tr><th id="L343"><a href="#L343">343</a></th><td></td></tr><tr><th id="L344"><a href="#L344">344</a></th><td>        if self.atoms:</td></tr><tr><th id="L345"><a href="#L345">345</a></th><td>            cell_atoms = self.atoms.get_cell()</td></tr><tr><th id="L346"><a href="#L346">346</a></th><td></td></tr><tr><th id="L347"><a href="#L347">347</a></th><td>            rotation_lammps2ase = np.dot(np.linalg.inv(np.array(cell)), cell_atoms)</td></tr><tr><th id="L348"><a href="#L348">348</a></th><td>#           print "rotation_lammps2ase:"</td></tr><tr><th id="L349"><a href="#L349">349</a></th><td>#           print rotation_lammps2ase</td></tr><tr><th id="L350"><a href="#L350">350</a></th><td>#           print np.transpose(rotation_lammps2ase)</td></tr><tr><th id="L351"><a href="#L351">351</a></th><td>#           print np.linalg.det(rotation_lammps2ase)                                            # check for orthogonality of matrix 'rotation'</td></tr><tr><th id="L352"><a href="#L352">352</a></th><td>#           print np.dot(rotation_lammps2ase, np.transpose(rotation_lammps2ase))                # check for orthogonality of matrix 'rotation'</td></tr><tr><th id="L353"><a href="#L353">353</a></th><td></td></tr><tr><th id="L354"><a href="#L354">354</a></th><td>            type_atoms = self.atoms.get_atomic_numbers()</td></tr><tr><th id="L355"><a href="#L355">355</a></th><td>            positions_atoms = np.array( [np.dot(np.array(r), rotation_lammps2ase) for r in positions] )</td></tr><tr><th id="L356"><a href="#L356">356</a></th><td>            velocities_atoms = np.array( [np.dot(np.array(v), rotation_lammps2ase) for v in velocities] )</td></tr><tr><th id="L357"><a href="#L357">357</a></th><td>            forces_atoms = np.array( [np.dot(np.array(f), rotation_lammps2ase) for f in forces] )</td></tr><tr><th id="L358"><a href="#L358">358</a></th><td></td></tr><tr><th id="L359"><a href="#L359">359</a></th><td>        if (set_atoms):</td></tr><tr><th id="L360"><a href="#L360">360</a></th><td>            # assume periodic boundary conditions here (like also below in write_lammps) &lt;- TODO:?!</td></tr><tr><th id="L361"><a href="#L361">361</a></th><td>            self.atoms = Atoms(type_atoms, positions=positions_atoms, cell=cell_atoms, pbc=True)</td></tr><tr><th id="L362"><a href="#L362">362</a></th><td></td></tr><tr><th id="L363"><a href="#L363">363</a></th><td>        self.forces = forces_atoms</td></tr><tr><th id="L364"><a href="#L364">364</a></th><td></td></tr><tr><th id="L365"><a href="#L365">365</a></th><td></td></tr><tr><th id="L366"><a href="#L366">366</a></th><td># could perhaps go into io.lammps in the future...</td></tr><tr><th id="L367"><a href="#L367">367</a></th><td>#from ase.parallel import paropen</td></tr><tr><th id="L368"><a href="#L368">368</a></th><td>def write_lammps(fileobj, atoms, force_skew=False):</td></tr><tr><th id="L369"><a href="#L369">369</a></th><td>    """Method which writes atomic structure data to a LAMMPS data file."""</td></tr><tr><th id="L370"><a href="#L370">370</a></th><td>    if isinstance(fileobj, file):</td></tr><tr><th id="L371"><a href="#L371">371</a></th><td>        f = fileobj</td></tr><tr><th id="L372"><a href="#L372">372</a></th><td>    elif isinstance(fileobj, str):</td></tr><tr><th id="L373"><a href="#L373">373</a></th><td>        f = paropen(fileobj, 'w')</td></tr><tr><th id="L374"><a href="#L374">374</a></th><td>    else :</td></tr><tr><th id="L375"><a href="#L375">375</a></th><td>        raise TypeError('fileobj is not of type file or str!')</td></tr><tr><th id="L376"><a href="#L376">376</a></th><td></td></tr><tr><th id="L377"><a href="#L377">377</a></th><td>    if isinstance(atoms, list):</td></tr><tr><th id="L378"><a href="#L378">378</a></th><td>        if len(atoms) &gt; 1:</td></tr><tr><th id="L379"><a href="#L379">379</a></th><td>            raise ValueError('Can only write one configuration to a lammps data file!')</td></tr><tr><th id="L380"><a href="#L380">380</a></th><td>        atoms = atoms[0]</td></tr><tr><th id="L381"><a href="#L381">381</a></th><td></td></tr><tr><th id="L382"><a href="#L382">382</a></th><td>    f.write(f.name + " (written by ASE) \n\n")</td></tr><tr><th id="L383"><a href="#L383">383</a></th><td></td></tr><tr><th id="L384"><a href="#L384">384</a></th><td>    symbols = atoms.get_chemical_symbols()</td></tr><tr><th id="L385"><a href="#L385">385</a></th><td>    n_atoms = len(symbols)</td></tr><tr><th id="L386"><a href="#L386">386</a></th><td>    f.write("%d \t atoms \n" % n_atoms)</td></tr><tr><th id="L387"><a href="#L387">387</a></th><td></td></tr><tr><th id="L388"><a href="#L388">388</a></th><td>    # Uniqify 'symbols' list and sort to a new list 'species'.</td></tr><tr><th id="L389"><a href="#L389">389</a></th><td>    # Sorting is important in case of multi-component systems:</td></tr><tr><th id="L390"><a href="#L390">390</a></th><td>    # This way it is assured that LAMMPS atom types are always</td></tr><tr><th id="L391"><a href="#L391">391</a></th><td>    # assigned predictively according to the alphabetic order </td></tr><tr><th id="L392"><a href="#L392">392</a></th><td>    # of the species present in the current system. </td></tr><tr><th id="L393"><a href="#L393">393</a></th><td>    # (Hence e.g. the mapping in interaction statements for alloy</td></tr><tr><th id="L394"><a href="#L394">394</a></th><td>    #  potentials depending on the order of the elements in the</td></tr><tr><th id="L395"><a href="#L395">395</a></th><td>    #  external potential can also be safely adjusted accordingly.)</td></tr><tr><th id="L396"><a href="#L396">396</a></th><td>    species = sorted(list(set(symbols)))</td></tr><tr><th id="L397"><a href="#L397">397</a></th><td>    n_atom_types = len(species)</td></tr><tr><th id="L398"><a href="#L398">398</a></th><td>    f.write("%d \t atom types \n" % n_atom_types)</td></tr><tr><th id="L399"><a href="#L399">399</a></th><td></td></tr><tr><th id="L400"><a href="#L400">400</a></th><td></td></tr><tr><th id="L401"><a href="#L401">401</a></th><td>    ### cell -&gt; bounding box for LAMMPS</td></tr><tr><th id="L402"><a href="#L402">402</a></th><td></td></tr><tr><th id="L403"><a href="#L403">403</a></th><td>    ## A) simple version</td></tr><tr><th id="L404"><a href="#L404">404</a></th><td>    # - For orthogonal cells which are properly aligned with the reference coordinate system only.</td></tr><tr><th id="L405"><a href="#L405">405</a></th><td>#    cell = atoms.get_cell()</td></tr><tr><th id="L406"><a href="#L406">406</a></th><td>#    lo = [0.0, 0.0, 0.0]</td></tr><tr><th id="L407"><a href="#L407">407</a></th><td>#    hi = sum(cell[0:3])</td></tr><tr><th id="L408"><a href="#L408">408</a></th><td>#    f.write("%15.8f %15.8f \t xlo xhi \n" % (lo[0], hi[0]))</td></tr><tr><th id="L409"><a href="#L409">409</a></th><td>#    f.write("%15.8f %15.8f \t ylo yhi \n" % (lo[1], hi[1]))</td></tr><tr><th id="L410"><a href="#L410">410</a></th><td>#    f.write("%15.8f %15.8f \t zlo zhi \n" % (lo[2], hi[2]))</td></tr><tr><th id="L411"><a href="#L411">411</a></th><td>#    f.write("\n\n")</td></tr><tr><th id="L412"><a href="#L412">412</a></th><td></td></tr><tr><th id="L413"><a href="#L413">413</a></th><td>    ## B) convert to (perhaps) triclinic cells in LAMMPS</td></tr><tr><th id="L414"><a href="#L414">414</a></th><td>    </td></tr><tr><th id="L415"><a href="#L415">415</a></th><td>    # B.1) calculate lengths of and angles between cell vectors</td></tr><tr><th id="L416"><a href="#L416">416</a></th><td>    #     (conventional crystallographic nomenclature)</td></tr><tr><th id="L417"><a href="#L417">417</a></th><td>    cell = atoms.get_cell()</td></tr><tr><th id="L418"><a href="#L418">418</a></th><td>    a = np.linalg.norm(cell[0])</td></tr><tr><th id="L419"><a href="#L419">419</a></th><td>    b = np.linalg.norm(cell[1])</td></tr><tr><th id="L420"><a href="#L420">420</a></th><td>    c = np.linalg.norm(cell[2])</td></tr><tr><th id="L421"><a href="#L421">421</a></th><td>    alpha = np.arccos( np.vdot(cell[1], cell[2]) / (b * c) )</td></tr><tr><th id="L422"><a href="#L422">422</a></th><td>    beta = np.arccos( np.vdot(cell[0], cell[2]) / (a * c) )</td></tr><tr><th id="L423"><a href="#L423">423</a></th><td>    gamma = np.arccos( np.vdot(cell[0], cell[1]) / (a * b) )</td></tr><tr><th id="L424"><a href="#L424">424</a></th><td>#    print [a, b, c]</td></tr><tr><th id="L425"><a href="#L425">425</a></th><td>#    print map(np.degrees, [alpha, beta, gamma])</td></tr><tr><th id="L426"><a href="#L426">426</a></th><td></td></tr><tr><th id="L427"><a href="#L427">427</a></th><td>    # B.2) construct bounding box edges and skew vector of </td></tr><tr><th id="L428"><a href="#L428">428</a></th><td>    #     corresponding (perhaps rotated!) LAMMPS cell</td></tr><tr><th id="L429"><a href="#L429">429</a></th><td>    # http://lammps.sandia.gov/doc/read_data.html :</td></tr><tr><th id="L430"><a href="#L430">430</a></th><td>    #   a_LAMMPS = (xhi-xlo,0,0); b_LAMMPS = (xy,yhi-ylo,0); c_LAMMPS = (xz,yz,zhi-zlo)</td></tr><tr><th id="L431"><a href="#L431">431</a></th><td>    xlo = ylo = zlo = 0.0</td></tr><tr><th id="L432"><a href="#L432">432</a></th><td>    # this choice of origin simplifies things:</td></tr><tr><th id="L433"><a href="#L433">433</a></th><td>    #   a_LAMMPS = (xhi,0,0); b_LAMMPS = (xy,yhi,0); c_LAMMPS = (xz,yz,zhi)</td></tr><tr><th id="L434"><a href="#L434">434</a></th><td>    xhi = a                     # a_LAMMPS</td></tr><tr><th id="L435"><a href="#L435">435</a></th><td>    xy = np.cos(gamma) * b      # b_LAMMPS</td></tr><tr><th id="L436"><a href="#L436">436</a></th><td>    yhi = np.sin(gamma) * b</td></tr><tr><th id="L437"><a href="#L437">437</a></th><td>    xz = np.cos(beta) * c       # c_LAMMPS</td></tr><tr><th id="L438"><a href="#L438">438</a></th><td>    yz = ( b * c * np.cos(alpha) - xy * xz ) / yhi</td></tr><tr><th id="L439"><a href="#L439">439</a></th><td>    zhi = np.sqrt( c**2 - xz**2 - yz**2 )</td></tr><tr><th id="L440"><a href="#L440">440</a></th><td>    </td></tr><tr><th id="L441"><a href="#L441">441</a></th><td>    # B.3) obtain rotation of original cell with respect to LAMMPS cell</td></tr><tr><th id="L442"><a href="#L442">442</a></th><td>    #     to properly tranform the atoms contained within (see below!)</td></tr><tr><th id="L443"><a href="#L443">443</a></th><td>    # IMPORTANT: use same convention as in cell of ASE atoms object, i.e.</td></tr><tr><th id="L444"><a href="#L444">444</a></th><td>    #            cell vectors are ROW VECTORS in numpy array</td></tr><tr><th id="L445"><a href="#L445">445</a></th><td>    cell_lammps = np.array([[xhi-xlo,0,0],[xy,yhi-ylo,0],[xz,yz,zhi-zlo]])</td></tr><tr><th id="L446"><a href="#L446">446</a></th><td>#    a_lammps = np.linalg.norm(cell_lammps[0])</td></tr><tr><th id="L447"><a href="#L447">447</a></th><td>#    b_lammps = np.linalg.norm(cell_lammps[1])</td></tr><tr><th id="L448"><a href="#L448">448</a></th><td>#    c_lammps = np.linalg.norm(cell_lammps[2])</td></tr><tr><th id="L449"><a href="#L449">449</a></th><td>#    alpha_lammps = np.arccos( np.vdot(cell_lammps[1], cell_lammps[2]) / (b_lammps * c_lammps) )</td></tr><tr><th id="L450"><a href="#L450">450</a></th><td>#    beta_lammps = np.arccos( np.vdot(cell_lammps[0], cell_lammps[2]) / (a_lammps * c_lammps) )</td></tr><tr><th id="L451"><a href="#L451">451</a></th><td>#    gamma_lammps = np.arccos( np.vdot(cell_lammps[0], cell_lammps[1]) / (a_lammps * b_lammps) )</td></tr><tr><th id="L452"><a href="#L452">452</a></th><td>#    print [a_lammps, b_lammps, c_lammps]</td></tr><tr><th id="L453"><a href="#L453">453</a></th><td>#    print map(np.degrees, [alpha_lammps, beta_lammps, gamma_lammps])</td></tr><tr><th id="L454"><a href="#L454">454</a></th><td>    # IMPORTANT: need vector-(rotation-)matrix product (instead of matrix-vector product) here,</td></tr><tr><th id="L455"><a href="#L455">455</a></th><td>    #            cell vectors are ROW VECTORS (also see above)</td></tr><tr><th id="L456"><a href="#L456">456</a></th><td>    rotation = np.dot(np.linalg.inv(cell), cell_lammps)</td></tr><tr><th id="L457"><a href="#L457">457</a></th><td>#    print rotation</td></tr><tr><th id="L458"><a href="#L458">458</a></th><td>#    print np.transpose(rotation)</td></tr><tr><th id="L459"><a href="#L459">459</a></th><td>#    print np.linalg.det(rotation)                              # check for orthogonality of matrix 'rotation'</td></tr><tr><th id="L460"><a href="#L460">460</a></th><td>#    print np.dot(rotation, np.transpose(rotation))             # check for orthogonality of matrix 'rotation'</td></tr><tr><th id="L461"><a href="#L461">461</a></th><td>#    print np.dot(rotation, cell[0])</td></tr><tr><th id="L462"><a href="#L462">462</a></th><td>#    print np.dot(rotation, cell[1])</td></tr><tr><th id="L463"><a href="#L463">463</a></th><td>#    print np.dot(rotation, cell[2])</td></tr><tr><th id="L464"><a href="#L464">464</a></th><td></td></tr><tr><th id="L465"><a href="#L465">465</a></th><td>    # B.4.1) write bounding box edges</td></tr><tr><th id="L466"><a href="#L466">466</a></th><td>    f.write("%15.8f %15.8f \t\t\t xlo xhi \n" % (xlo, xhi))</td></tr><tr><th id="L467"><a href="#L467">467</a></th><td>    f.write("%15.8f %15.8f \t\t\t ylo yhi \n" % (ylo, yhi))</td></tr><tr><th id="L468"><a href="#L468">468</a></th><td>    f.write("%15.8f %15.8f \t\t\t zlo zhi \n" % (zlo, zhi))</td></tr><tr><th id="L469"><a href="#L469">469</a></th><td>    </td></tr><tr><th id="L470"><a href="#L470">470</a></th><td>    # B.4.2) sanitize and write skew vector (if necessary)</td></tr><tr><th id="L471"><a href="#L471">471</a></th><td>    # The too simple version </td></tr><tr><th id="L472"><a href="#L472">472</a></th><td>    #    f.write("%20.10f %20.10f %20.10f \t xy xz yz \n" % (xy, xz, yz))    </td></tr><tr><th id="L473"><a href="#L473">473</a></th><td>    # can make LAMMPS (easily) reject the definition of a triclinic box because</td></tr><tr><th id="L474"><a href="#L474">474</a></th><td>    #   a) skew vector components calculated above are outside the respective ranges </td></tr><tr><th id="L475"><a href="#L475">475</a></th><td>    #      which LAMMPS expects as described here</td></tr><tr><th id="L476"><a href="#L476">476</a></th><td>    #           http://lammps.sandia.gov/doc/read_data.html</td></tr><tr><th id="L477"><a href="#L477">477</a></th><td>    #      and coded here</td></tr><tr><th id="L478"><a href="#L478">478</a></th><td>    #           domain.cpp -&gt; Domain::set_initial_box()</td></tr><tr><th id="L479"><a href="#L479">479</a></th><td>    #   b) rounding up in the skew vector output can make the result violate those </td></tr><tr><th id="L480"><a href="#L480">480</a></th><td>    #      ranges for purely numeric reasons (perhaps even more frequent problem than a)</td></tr><tr><th id="L481"><a href="#L481">481</a></th><td>    #      and definitely more stupid!)</td></tr><tr><th id="L482"><a href="#L482">482</a></th><td>    # More elaborate solution addressing the issues described above:</td></tr><tr><th id="L483"><a href="#L483">483</a></th><td>    # -&gt; a):</td></tr><tr><th id="L484"><a href="#L484">484</a></th><td>    # Fold back skew vector components calculated above to equivalent ones complying with</td></tr><tr><th id="L485"><a href="#L485">485</a></th><td>    # what LAMMPS expects (see references above) and describing another ("less skewed") </td></tr><tr><th id="L486"><a href="#L486">486</a></th><td>    # unit cell for the same lattice:</td></tr><tr><th id="L487"><a href="#L487">487</a></th><td>    #           t -&gt; t_folded in [-period/2, +period/2]</td></tr><tr><th id="L488"><a href="#L488">488</a></th><td>    def fold_skew_dec(t, period) :</td></tr><tr><th id="L489"><a href="#L489">489</a></th><td>        t_dec = dec.Decimal(repr(t))</td></tr><tr><th id="L490"><a href="#L490">490</a></th><td>        period_dec = dec.Decimal(repr(period))</td></tr><tr><th id="L491"><a href="#L491">491</a></th><td>        # dec.ROUND_HALF_DOWN leaves all t in [-period/2, +period/2] unchanged whereas</td></tr><tr><th id="L492"><a href="#L492">492</a></th><td>        # dec.ROUND_HALF_UP does the same for all t in (-period/2, +period/2) but</td></tr><tr><th id="L493"><a href="#L493">493</a></th><td>        # swaps the boundary values: </td></tr><tr><th id="L494"><a href="#L494">494</a></th><td>        #       t=-period/2 -&gt; t_folded=+period/2</td></tr><tr><th id="L495"><a href="#L495">495</a></th><td>        #       t=+period/2 -&gt; t_folded=-period/2</td></tr><tr><th id="L496"><a href="#L496">496</a></th><td>        t_folded_dec = t_dec - period_dec * (t_dec/period_dec).to_integral_value(dec.ROUND_HALF_DOWN)</td></tr><tr><th id="L497"><a href="#L497">497</a></th><td>        return t_folded_dec</td></tr><tr><th id="L498"><a href="#L498">498</a></th><td>    skew_folded_dec = [ fold_skew_dec(xy, xhi-xlo), fold_skew_dec(xz, xhi-xlo), fold_skew_dec(yz, yhi-ylo) ]</td></tr><tr><th id="L499"><a href="#L499">499</a></th><td>    # -&gt; b):</td></tr><tr><th id="L500"><a href="#L500">500</a></th><td>    # This should be kept in sync with the accuracy used for writing the bounding box edges above</td></tr><tr><th id="L501"><a href="#L501">501</a></th><td>    prec_dec = dec.Decimal('1E-8')</td></tr><tr><th id="L502"><a href="#L502">502</a></th><td>    # Make sure to always round down and hence avoid numerical problems with skew vector in LAMMPS.</td></tr><tr><th id="L503"><a href="#L503">503</a></th><td>    skew_folded_and_rounded_dec = [ d.quantize(prec_dec, dec.ROUND_DOWN) for d in skew_folded_dec ]</td></tr><tr><th id="L504"><a href="#L504">504</a></th><td>    # Do not write zero skew vector. </td></tr><tr><th id="L505"><a href="#L505">505</a></th><td>    # (I.e. calculate orthogonal cell with LAMMPS in that case!)</td></tr><tr><th id="L506"><a href="#L506">506</a></th><td>    has_skew = any( [not d.is_zero() for d in skew_folded_and_rounded_dec] )</td></tr><tr><th id="L507"><a href="#L507">507</a></th><td>    if force_skew or has_skew:</td></tr><tr><th id="L508"><a href="#L508">508</a></th><td>        f.write( "%15s %15s %15s \t xy xz yz \n" % tuple( str(d) for d in skew_folded_and_rounded_dec ) )</td></tr><tr><th id="L509"><a href="#L509">509</a></th><td>    f.write("\n\n")</td></tr><tr><th id="L510"><a href="#L510">510</a></th><td></td></tr><tr><th id="L511"><a href="#L511">511</a></th><td></td></tr><tr><th id="L512"><a href="#L512">512</a></th><td>    ### atoms</td></tr><tr><th id="L513"><a href="#L513">513</a></th><td></td></tr><tr><th id="L514"><a href="#L514">514</a></th><td>    ## A) simple version</td></tr><tr><th id="L515"><a href="#L515">515</a></th><td>    # - For orthogonal cells which are properly aligned with the reference coordinate system only.</td></tr><tr><th id="L516"><a href="#L516">516</a></th><td>#    f.write("Atoms \n\n")</td></tr><tr><th id="L517"><a href="#L517">517</a></th><td>#    for i, (x, y, z) in enumerate(atoms.get_positions()):</td></tr><tr><th id="L518"><a href="#L518">518</a></th><td>#        s = species.index(symbols[i]) + 1</td></tr><tr><th id="L519"><a href="#L519">519</a></th><td>#        f.write("%6d %3d %22.15f %22.15f %22.15f \n" % (i+1, s, x, y, z))</td></tr><tr><th id="L520"><a href="#L520">520</a></th><td></td></tr><tr><th id="L521"><a href="#L521">521</a></th><td>    ## B) convert to (perhaps) triclinic cells in LAMMPS:</td></tr><tr><th id="L522"><a href="#L522">522</a></th><td>    # - Apply rotation determined above for this case.</td></tr><tr><th id="L523"><a href="#L523">523</a></th><td>    # - Lattice translation in case of a folded skew vector above should not be necessary,</td></tr><tr><th id="L524"><a href="#L524">524</a></th><td>    #   since the resulting (perhaps) new unit cell does describe the same lattice.</td></tr><tr><th id="L525"><a href="#L525">525</a></th><td>    #   Therefore the (rotated) atoms are still at correct lattice sites, only some of them</td></tr><tr><th id="L526"><a href="#L526">526</a></th><td>    #   probably are outside of that new unit cell, which LAMMPS can hopefully deal with.</td></tr><tr><th id="L527"><a href="#L527">527</a></th><td>    f.write("Atoms \n\n")</td></tr><tr><th id="L528"><a href="#L528">528</a></th><td>    for i, r in enumerate(atoms.get_positions()):</td></tr><tr><th id="L529"><a href="#L529">529</a></th><td>        s = species.index(symbols[i]) + 1</td></tr><tr><th id="L530"><a href="#L530">530</a></th><td>        [x, y, z] = np.dot(r, rotation)</td></tr><tr><th id="L531"><a href="#L531">531</a></th><td>#        print r, [x, y, z]</td></tr><tr><th id="L532"><a href="#L532">532</a></th><td>        f.write("%6d %3d %22.15f %22.15f %22.15f \n" % (i+1, s, x, y, z))</td></tr><tr><th id="L533"><a href="#L533">533</a></th><td></td></tr><tr><th id="L534"><a href="#L534">534</a></th><td>    f.close()</td></tr><tr><th id="L535"><a href="#L535">535</a></th><td></td></tr><tr><th id="L536"><a href="#L536">536</a></th><td></td></tr><tr><th id="L537"><a href="#L537">537</a></th><td>if __name__ == "__main__":</td></tr><tr><th id="L538"><a href="#L538">538</a></th><td></td></tr><tr><th id="L539"><a href="#L539">539</a></th><td>    pair_style = "eam"</td></tr><tr><th id="L540"><a href="#L540">540</a></th><td>    Pd_eam_file = "Pd_u3.eam"</td></tr><tr><th id="L541"><a href="#L541">541</a></th><td>    pair_coeff = [ "* * " + Pd_eam_file ]</td></tr><tr><th id="L542"><a href="#L542">542</a></th><td>    parameters = { "pair_style" : pair_style, "pair_coeff" : pair_coeff }</td></tr><tr><th id="L543"><a href="#L543">543</a></th><td>    files = [ Pd_eam_file ]</td></tr><tr><th id="L544"><a href="#L544">544</a></th><td>    calc = LAMMPS(parameters=parameters, files=files)</td></tr><tr><th id="L545"><a href="#L545">545</a></th><td></td></tr><tr><th id="L546"><a href="#L546">546</a></th><td>    from ase import Atoms</td></tr><tr><th id="L547"><a href="#L547">547</a></th><td>    a0 = 3.93</td></tr><tr><th id="L548"><a href="#L548">548</a></th><td>    b0 = a0 / 2.0</td></tr><tr><th id="L549"><a href="#L549">549</a></th><td>    if True:</td></tr><tr><th id="L550"><a href="#L550">550</a></th><td>        bulk = Atoms(['Pd']*4,</td></tr><tr><th id="L551"><a href="#L551">551</a></th><td>                     positions=[(0,0,0),(b0,b0,0),(b0,0,b0),(0,b0,b0)],</td></tr><tr><th id="L552"><a href="#L552">552</a></th><td>                     cell=[a0]*3,</td></tr><tr><th id="L553"><a href="#L553">553</a></th><td>                     pbc=True)</td></tr><tr><th id="L554"><a href="#L554">554</a></th><td>        # test get_forces</td></tr><tr><th id="L555"><a href="#L555">555</a></th><td>        print "forces for a = %f" % a0</td></tr><tr><th id="L556"><a href="#L556">556</a></th><td>        print calc.get_forces(bulk)</td></tr><tr><th id="L557"><a href="#L557">557</a></th><td>        # single points for various lattice constants</td></tr><tr><th id="L558"><a href="#L558">558</a></th><td>        bulk.set_calculator(calc)</td></tr><tr><th id="L559"><a href="#L559">559</a></th><td>        for n in range(-5,5,1):</td></tr><tr><th id="L560"><a href="#L560">560</a></th><td>            a = a0 * (1 + n/100.0)</td></tr><tr><th id="L561"><a href="#L561">561</a></th><td>            bulk.set_cell([a]*3)</td></tr><tr><th id="L562"><a href="#L562">562</a></th><td>            print "a : %f , total energy : %f" % (a, bulk.get_potential_energy())</td></tr><tr><th id="L563"><a href="#L563">563</a></th><td></td></tr><tr><th id="L564"><a href="#L564">564</a></th><td>#    calc.clean(force=True)</td></tr><tr><th id="L565"><a href="#L565">565</a></th><td>    calc.clean()</td></tr><tr><th id="L566"><a href="#L566">566</a></th><td></td></tr><tr><th id="L567"><a href="#L567">567</a></th><td></td></tr><tr><th id="L568"><a href="#L568">568</a></th><td></td></tr></tbody></table>
+        for f in files:
+            shutil.copy(f, os.path.join(self.dir, f))
 
-      </div>
-      <div id="help"><strong>Note:</strong> See <a href="/projects/ase/wiki/TracBrowser">TracBrowser</a>
-        for help on using the repository browser.</div>
-      <div id="anydiff">
-        <form action="/projects/ase/diff" method="get">
-          <div class="buttons">
-            <input type="hidden" name="new_path" value="/trunk/ase/calculators/lammps.py" />
-            <input type="hidden" name="old_path" value="/trunk/ase/calculators/lammps.py" />
-            <input type="hidden" name="new_rev" />
-            <input type="hidden" name="old_rev" />
-            <input type="submit" value="View changes..." title="Select paths and revs for Diff" />
-          </div>
-        </form>
-      </div>
-    </div>
-    <div id="altlinks">
-      <h3>Download in other formats:</h3>
-      <ul>
-        <li class="first">
-          <a rel="nofollow" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?format=txt">Plain Text</a>
-        </li><li class="last">
-          <a rel="nofollow" href="/projects/ase/export/2073/trunk/ase/calculators/lammps.py">Original Format</a>
-        </li>
-      </ul>
-    </div>
-    </div>
-    <div id="footer" lang="en" xml:lang="en"><hr />
-      <a id="tracpowered" href="http://trac.edgewall.org/"><img src="/projects/ase/chrome/common/trac_logo_mini.png" height="30" width="107" alt="Trac Powered" /></a>
-      <p class="left">Powered by <a href="/projects/ase/about"><strong>Trac 0.12</strong></a><br />
-        By <a href="http://www.edgewall.org/">Edgewall Software</a>.</p>
-      <p class="right">Visit the Trac open source project at<br /><a href="http://trac.edgewall.org/">http://trac.edgewall.org/</a></p>
-    </div>
-  </body>
-</html>
\ No newline at end of file
+    def clean(self, force=False):
+        if os.path.isdir(self.dir):
+            files_in_dir = os.listdir(self.dir)
+            files_copied = self.files
+            if (len(files_in_dir) == 0):
+                os.rmdir(self.dir)
+            elif (set(files_in_dir).issubset(set(files_copied)) or force):
+                shutil.rmtree(self.dir)
+
+    def get_potential_energy(self, atoms):
+        self.update(atoms)
+        return self.potential_energy
+
+    def get_forces(self, atoms):
+        self.update(atoms)
+        return self.forces.copy()
+
+    def get_stress(self, atoms):
+        raise NotImplementedError
+
+    def update(self, atoms):
+        # TODO: check if (re-)calculation is necessary
+        self.calculate(atoms)
+
+    def calculate(self, atoms):
+        self.atoms = atoms.copy()
+        self.run()
+
+    def run(self):
+        """Method which explicitely runs LAMMPS."""
+
+        self.calls += 1
+
+        # set LAMMPS command from environment variable
+        if os.environ.has_key('LAMMPS_COMMAND'):
+            lammps_command = os.path.abspath(os.environ['LAMMPS_COMMAND'])
+        else:
+            self.clean()
+            raise RuntimeError('Please set LAMMPS_COMMAND environment variable')
+        if os.environ.has_key('LAMMPS_OPTIONS'):
+            lammps_options = os.environ['LAMMPS_OPTIONS']
+        else:
+            lammps_options = "-echo log -screen none"
+
+        # change into subdirectory for LAMMPS calculations
+        cwd = os.getcwd()
+        os.chdir(self.dir)
+
+        # setup file names for LAMMPS calculation
+        label = "%s-%06d" % (self.label, self.calls)
+        lammps_data = "data." + label
+        lammps_in = "in." + label
+        lammps_trj = label + ".lammpstrj"
+        lammps_log = label + ".log"
+
+        # write LAMMPS input files
+        self.write_lammps_data(lammps_data=lammps_data)
+        self.write_lammps_in(lammps_in=lammps_in, lammps_data=lammps_data, lammps_trj=lammps_trj, parameters=self.parameters)
+
+        # run LAMMPS
+        # TODO: check for successful completion (based on output files?!) of LAMMPS call...
+        exitcode = os.system("%s %s -in %s -log %s" % (lammps_command, lammps_options, lammps_in, lammps_log))
+        if exitcode != 0:
+            cwd = os.getcwd()
+            raise RuntimeError("LAMMPS exited in %s with exit code: %d.  " % (cwd,exitcode))
+
+        # read LAMMPS output files
+        self.read_lammps_log(lammps_log=lammps_log)
+        self.read_lammps_trj(lammps_trj=lammps_trj)
+
+        # change back to previous working directory
+        os.chdir(cwd)
+
+    def write_lammps_data(self, lammps_data=None):
+        """Method which writes a LAMMPS data file with atomic structure."""
+        if (lammps_data == None):
+            lammps_data = "data." + self.label
+        write_lammps(lammps_data, self.atoms, force_skew=self.always_triclinic)
+
+    def write_lammps_in(self, lammps_in=None, lammps_data=None, lammps_trj=None, parameters={}):
+        """Method which writes a LAMMPS in file with run parameters and settings."""
+        if (lammps_in == None):
+            lammps_in = "in." + self.label
+        if (lammps_data == None):
+            lammps_data = "data." + self.label
+        if (lammps_trj == None):
+            lammps_trj = self.label + ".lammpstrj"
+
+        f = paropen(lammps_in, 'w')
+        f.write("# " + f.name + " (written by ASE) \n")
+        f.write("\n")
+
+        f.write("### variables \n")
+        f.write("variable data_file index \"%s\" \n" % lammps_data)
+        f.write("variable dump_file index \"%s\" \n" % lammps_trj)
+        f.write("\n\n")
+
+        pbc = self.atoms.get_pbc()
+        f.write("### simulation box \n")
+        f.write("units metal \n")
+        if ("boundary" in parameters):
+            f.write("boundary %s \n" % parameters["boundary"])
+        else:
+            f.write("boundary %c %c %c \n" % tuple('sp'[x] for x in pbc))
+        f.write("atom_modify sort 0 0.0 \n")
+        for key in ("neighbor" ,"newton"):
+            if key in parameters:
+                f.write("%s %s \n" % (key, parameters[key]))
+        f.write("\n")
+        f.write("read_data ${data_file} \n")
+        f.write("\n\n")
+
+        f.write("### interactions \n")
+        if ( ("pair_style" in parameters) and ("pair_coeff" in parameters) ):
+            pair_style = parameters["pair_style"]
+            f.write("pair_style %s \n" % pair_style)
+            for pair_coeff in parameters["pair_coeff"]:
+                f.write("pair_coeff %s \n" % pair_coeff)
+            if "mass" in parameters:
+                for mass in parameters["mass"]:
+                    f.write("mass %s \n" % mass)
+        else:
+            # simple default parameters that should always make the LAMMPS calculation run
+            #    pair_style 	lj/cut 2.5
+            #    pair_coeff 	* * 1 1
+            #    mass           * 1.0        else:
+            f.write("pair_style lj/cut 2.5 \n")
+            f.write("pair_coeff * * 1 1 \n")
+            f.write("mass * 1.0 \n")
+        f.write("\n")
+
+        f.write("### run \n")
+        f.write("fix fix_nve all nve \n")
+        f.write("\n")
+        f.write("dump dump_all all custom 1 ${dump_file} id type x y z vx vy vz fx fy fz \n")
+        f.write("\n")
+        f.write("thermo_style custom step temp ke pe etotal cpu \n")
+        f.write("thermo_modify format 1 %4d format 2 %9.3f format 3 %20.6f format 4 %20.6f format 5 %20.6f format 6 %9.3f \n")
+        f.write("thermo 1 \n")
+        f.write("\n")
+        if ("minimize" in parameters):
+            f.write("minimize %s \n" % parameters["minimize"])
+        if ("run" in parameters):
+            f.write("run %s \n" % parameters["run"])
+        if not ( ("minimize" in parameters) or ("run" in parameters) ):
+            f.write("run 0 \n")
+
+        f.close()
+
+    def read_lammps_log(self, lammps_log=None):
+        """Method which reads a LAMMPS output log file."""
+        if (lammps_log == None):
+            lammps_log = self.label + ".log"
+
+        epot = 0.0
+
+        f = paropen(lammps_log, 'r')
+        while True:
+            line = f.readline()
+            if not line:
+                break
+            # get potential energy of first step (if more are done)
+            if "PotEng" in line:
+                i = line.split().index("PotEng")
+                line = f.readline()
+                epot = float(line.split()[i])
+        f.close()
+
+#        print epot
+
+        self.potential_energy = epot
+
+    def read_lammps_trj(self, lammps_trj=None, set_atoms=False):
+        """Method which reads a LAMMPS dump file."""
+        if (lammps_trj == None):
+            lammps_trj = self.label + ".lammpstrj"
+
+        f = paropen(lammps_trj, 'r')
+        while True:
+            line = f.readline()
+
+            if not line:
+                break
+
+            #TODO: extend to proper dealing with multiple steps in one trajectory file
+            if "ITEM: TIMESTEP" in line:
+                n_atoms = 0
+                lo = [] ; hi = [] ; tilt = []
+                id = [] ; type = []
+                positions = [] ; velocities = [] ; forces = []
+
+            if "ITEM: NUMBER OF ATOMS" in line:
+                line = f.readline()
+                n_atoms = int(line.split()[0])
+            
+            if "ITEM: BOX BOUNDS" in line:
+                # save labels behind "ITEM: BOX BOUNDS" in triclinic case (>=lammps-7Jul09)
+                tilt_items = line.split()[3:]
+                for i in range(3):
+                    line = f.readline()
+                    fields = line.split()
+                    lo.append(float(fields[0]))
+                    hi.append(float(fields[1]))
+                    if (len(fields) >= 3):
+                        tilt.append(float(fields[2]))
+            
+            if "ITEM: ATOMS" in line:
+                # (reliably) identify values by labels behind "ITEM: ATOMS" - requires >=lammps-7Jul09
+                # create corresponding index dictionary before iterating over atoms to (hopefully) speed up lookups...
+                atom_attributes = {}
+                for (i, x) in enumerate(line.split()[2:]):
+                    atom_attributes[x] = i
+                for n in range(n_atoms):
+                    line = f.readline()
+                    fields = line.split()
+                    id.append( int(fields[atom_attributes['id']]) )
+                    type.append( int(fields[atom_attributes['type']]) )
+                    positions.append( [ float(fields[atom_attributes[x]]) for x in ['x', 'y', 'z'] ] )
+                    velocities.append( [ float(fields[atom_attributes[x]]) for x in ['vx', 'vy', 'vz'] ] )
+                    forces.append( [ float(fields[atom_attributes[x]]) for x in ['fx', 'fy', 'fz'] ] )
+        f.close()
+
+        # determine cell tilt (triclinic case!)
+        if (len(tilt) >= 3):
+            # for >=lammps-7Jul09 use labels behind "ITEM: BOX BOUNDS" to assign tilt (vector) elements ...
+            if (len(tilt_items) >= 3):
+                xy = tilt[tilt_items.index('xy')]
+                xz = tilt[tilt_items.index('xz')]
+                yz = tilt[tilt_items.index('yz')]
+            # ... otherwise assume default order in 3rd column (if the latter was present)
+            else:
+                xy = tilt[0]
+                xz = tilt[1]
+                yz = tilt[2]
+        else:
+            xy = xz = yz = 0
+        xhilo = (hi[0] - lo[0]) - xy - xz
+        yhilo = (hi[1] - lo[1]) - yz
+        zhilo = (hi[2] - lo[2])
+	
+	# The simulation box bounds are included in each snapshot and if the box is triclinic (non-orthogonal), 
+	# then the tilt factors are also printed; see the region prism command for a description of tilt factors. 
+	# For triclinic boxes the box bounds themselves (first 2 quantities on each line) are a true "bounding box" 
+	# around the simulation domain, which means they include the effect of any tilt.
+	# [ http://lammps.sandia.gov/doc/dump.html , lammps-7Jul09 ]
+	#
+	# This *should* extract the lattice vectors that LAMMPS uses from the true "bounding box" printed in the dump file
+	# It might fail in some cases (negative tilts?!) due to the MIN / MAX construction of these box corners:
+	#
+	#	void Domain::set_global_box() 
+	#	[...]
+	#	  if (triclinic) {
+	#	    [...]
+	#	    boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy);
+	#	    boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz);
+	#	    boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz);
+	#	    boxlo_bound[2] = boxlo[2];
+	#
+	#	    boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy);
+	#	    boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz);
+	#	    boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);
+	#	    boxhi_bound[2] = boxhi[2];
+	#	  }
+	# [ lammps-7Jul09/src/domain.cpp ]
+	#
+        cell = [[xhilo,0,0],[xy,yhilo,0],[xz,yz,zhilo]]
+
+#       print n_atoms
+#       print lo
+#       print hi
+#       print id
+#       print type
+#       print positions
+#       print velocities
+#       print forces
+
+        # assume that LAMMPS does not reorder atoms internally
+#        print np.array(positions)
+#        print self.atoms.get_positions()
+#        print np.array(positions) - self.atoms.get_positions()
+
+        cell_atoms = np.array(cell)
+        type_atoms = np.array(type)
+#      velocities_atoms = np.array(velocities)
+        positions_atoms = np.array(positions)
+        forces_atoms = np.array(forces)
+
+        if self.atoms:
+            cell_atoms = self.atoms.get_cell()
+
+            rotation_lammps2ase = np.dot(np.linalg.inv(np.array(cell)), cell_atoms)
+#           print "rotation_lammps2ase:"
+#           print rotation_lammps2ase
+#           print np.transpose(rotation_lammps2ase)
+#           print np.linalg.det(rotation_lammps2ase)                                            # check for orthogonality of matrix 'rotation'
+#           print np.dot(rotation_lammps2ase, np.transpose(rotation_lammps2ase))                # check for orthogonality of matrix 'rotation'
+
+            type_atoms = self.atoms.get_atomic_numbers()
+            positions_atoms = np.array( [np.dot(np.array(r), rotation_lammps2ase) for r in positions] )
+            velocities_atoms = np.array( [np.dot(np.array(v), rotation_lammps2ase) for v in velocities] )
+            forces_atoms = np.array( [np.dot(np.array(f), rotation_lammps2ase) for f in forces] )
+
+        if (set_atoms):
+            # assume periodic boundary conditions here (like also below in write_lammps) <- TODO:?!
+            self.atoms = Atoms(type_atoms, positions=positions_atoms, cell=cell_atoms, pbc=True)
+
+        self.forces = forces_atoms
+
+
+# could perhaps go into io.lammps in the future...
+#from ase.parallel import paropen
+def write_lammps(fileobj, atoms, force_skew=False):
+    """Method which writes atomic structure data to a LAMMPS data file."""
+    if isinstance(fileobj, file):
+        f = fileobj
+    elif isinstance(fileobj, str):
+        f = paropen(fileobj, 'w')
+    else :
+        raise TypeError('fileobj is not of type file or str!')
+
+    if isinstance(atoms, list):
+        if len(atoms) > 1:
+            raise ValueError('Can only write one configuration to a lammps data file!')
+        atoms = atoms[0]
+
+    f.write(f.name + " (written by ASE) \n\n")
+
+    symbols = atoms.get_chemical_symbols()
+    n_atoms = len(symbols)
+    f.write("%d \t atoms \n" % n_atoms)
+
+    # Uniqify 'symbols' list and sort to a new list 'species'.
+    # Sorting is important in case of multi-component systems:
+    # This way it is assured that LAMMPS atom types are always
+    # assigned predictively according to the alphabetic order 
+    # of the species present in the current system. 
+    # (Hence e.g. the mapping in interaction statements for alloy
+    #  potentials depending on the order of the elements in the
+    #  external potential can also be safely adjusted accordingly.)
+    species = sorted(list(set(symbols)))
+    n_atom_types = len(species)
+    f.write("%d \t atom types \n" % n_atom_types)
+
+
+    ### cell -> bounding box for LAMMPS
+
+    ## A) simple version
+    # - For orthogonal cells which are properly aligned with the reference coordinate system only.
+#    cell = atoms.get_cell()
+#    lo = [0.0, 0.0, 0.0]
+#    hi = sum(cell[0:3])
+#    f.write("%15.8f %15.8f \t xlo xhi \n" % (lo[0], hi[0]))
+#    f.write("%15.8f %15.8f \t ylo yhi \n" % (lo[1], hi[1]))
+#    f.write("%15.8f %15.8f \t zlo zhi \n" % (lo[2], hi[2]))
+#    f.write("\n\n")
+
+    ## B) convert to (perhaps) triclinic cells in LAMMPS
+    
+    # B.1) calculate lengths of and angles between cell vectors
+    #     (conventional crystallographic nomenclature)
+    cell = atoms.get_cell()
+    a = np.linalg.norm(cell[0])
+    b = np.linalg.norm(cell[1])
+    c = np.linalg.norm(cell[2])
+    alpha = np.arccos( np.vdot(cell[1], cell[2]) / (b * c) )
+    beta = np.arccos( np.vdot(cell[0], cell[2]) / (a * c) )
+    gamma = np.arccos( np.vdot(cell[0], cell[1]) / (a * b) )
+#    print [a, b, c]
+#    print map(np.degrees, [alpha, beta, gamma])
+
+    # B.2) construct bounding box edges and skew vector of 
+    #     corresponding (perhaps rotated!) LAMMPS cell
+    # http://lammps.sandia.gov/doc/read_data.html :
+    # 	a_LAMMPS = (xhi-xlo,0,0); b_LAMMPS = (xy,yhi-ylo,0); c_LAMMPS = (xz,yz,zhi-zlo)
+    xlo = ylo = zlo = 0.0
+    # this choice of origin simplifies things:
+    # 	a_LAMMPS = (xhi,0,0); b_LAMMPS = (xy,yhi,0); c_LAMMPS = (xz,yz,zhi)
+    xhi = a			# a_LAMMPS
+    xy = np.cos(gamma) * b	# b_LAMMPS
+    yhi = np.sin(gamma) * b
+    xz = np.cos(beta) * c	# c_LAMMPS
+    yz = ( b * c * np.cos(alpha) - xy * xz ) / yhi
+    zhi = np.sqrt( c**2 - xz**2 - yz**2 )
+    
+    # B.3) obtain rotation of original cell with respect to LAMMPS cell
+    #     to properly tranform the atoms contained within (see below!)
+    # IMPORTANT: use same convention as in cell of ASE atoms object, i.e.
+    #            cell vectors are ROW VECTORS in numpy array
+    cell_lammps = np.array([[xhi-xlo,0,0],[xy,yhi-ylo,0],[xz,yz,zhi-zlo]])
+#    a_lammps = np.linalg.norm(cell_lammps[0])
+#    b_lammps = np.linalg.norm(cell_lammps[1])
+#    c_lammps = np.linalg.norm(cell_lammps[2])
+#    alpha_lammps = np.arccos( np.vdot(cell_lammps[1], cell_lammps[2]) / (b_lammps * c_lammps) )
+#    beta_lammps = np.arccos( np.vdot(cell_lammps[0], cell_lammps[2]) / (a_lammps * c_lammps) )
+#    gamma_lammps = np.arccos( np.vdot(cell_lammps[0], cell_lammps[1]) / (a_lammps * b_lammps) )
+#    print [a_lammps, b_lammps, c_lammps]
+#    print map(np.degrees, [alpha_lammps, beta_lammps, gamma_lammps])
+    # IMPORTANT: need vector-(rotation-)matrix product (instead of matrix-vector product) here,
+    #            cell vectors are ROW VECTORS (also see above)
+    rotation = np.dot(np.linalg.inv(cell), cell_lammps)
+#    print rotation
+#    print np.transpose(rotation)
+#    print np.linalg.det(rotation)				# check for orthogonality of matrix 'rotation'
+#    print np.dot(rotation, np.transpose(rotation))		# check for orthogonality of matrix 'rotation'
+#    print np.dot(rotation, cell[0])
+#    print np.dot(rotation, cell[1])
+#    print np.dot(rotation, cell[2])
+
+    # B.4.1) write bounding box edges
+    f.write("%15.8f %15.8f \t\t\t xlo xhi \n" % (xlo, xhi))
+    f.write("%15.8f %15.8f \t\t\t ylo yhi \n" % (ylo, yhi))
+    f.write("%15.8f %15.8f \t\t\t zlo zhi \n" % (zlo, zhi))
+    
+    # B.4.2) sanitize and write skew vector (if necessary)
+    # The too simple version 
+    #    f.write("%20.10f %20.10f %20.10f \t xy xz yz \n" % (xy, xz, yz))    
+    # can make LAMMPS (easily) reject the definition of a triclinic box because
+    #	a) skew vector components calculated above are outside the respective ranges 
+    #	   which LAMMPS expects as described here
+    #		http://lammps.sandia.gov/doc/read_data.html
+    #	   and coded here
+    #		domain.cpp -> Domain::set_initial_box()
+    # 	b) rounding up in the skew vector output can make the result violate those 
+    #      ranges for purely numeric reasons (perhaps even more frequent problem than a)
+    #      and definitely more stupid!)
+    # More elaborate solution addressing the issues described above:
+    # -> a):
+    # Fold back skew vector components calculated above to equivalent ones complying with
+    # what LAMMPS expects (see references above) and describing another ("less skewed") 
+    # unit cell for the same lattice:
+    #		t -> t_folded in [-period/2, +period/2]
+    def fold_skew_dec(t, period) :
+        t_dec = dec.Decimal(repr(t))
+        period_dec = dec.Decimal(repr(period))
+	# dec.ROUND_HALF_DOWN leaves all t in [-period/2, +period/2] unchanged whereas
+	# dec.ROUND_HALF_UP does the same for all t in (-period/2, +period/2) but
+	# swaps the boundary values: 
+	#	t=-period/2 -> t_folded=+period/2
+	#	t=+period/2 -> t_folded=-period/2
+        t_folded_dec = t_dec - period_dec * (t_dec/period_dec).to_integral_value(dec.ROUND_HALF_DOWN)
+        return t_folded_dec
+    skew_folded_dec = [ fold_skew_dec(xy, xhi-xlo), fold_skew_dec(xz, xhi-xlo), fold_skew_dec(yz, yhi-ylo) ]
+    # -> b):
+    # This should be kept in sync with the accuracy used for writing the bounding box edges above
+    prec_dec = dec.Decimal('1E-8')
+    # Make sure to always round down and hence avoid numerical problems with skew vector in LAMMPS.
+    skew_folded_and_rounded_dec = [ d.quantize(prec_dec, dec.ROUND_DOWN) for d in skew_folded_dec ]
+    # Do not write zero skew vector. 
+    # (I.e. calculate orthogonal cell with LAMMPS in that case!)
+    has_skew = any( [not d.is_zero() for d in skew_folded_and_rounded_dec] )
+    if force_skew or has_skew:
+        f.write( "%15s %15s %15s \t xy xz yz \n" % tuple( str(d) for d in skew_folded_and_rounded_dec ) )
+    f.write("\n\n")
+
+
+    ### atoms
+
+    ## A) simple version
+    # - For orthogonal cells which are properly aligned with the reference coordinate system only.
+#    f.write("Atoms \n\n")
+#    for i, (x, y, z) in enumerate(atoms.get_positions()):
+#        s = species.index(symbols[i]) + 1
+#        f.write("%6d %3d %22.15f %22.15f %22.15f \n" % (i+1, s, x, y, z))
+
+    ## B) convert to (perhaps) triclinic cells in LAMMPS:
+    # - Apply rotation determined above for this case.
+    # - Lattice translation in case of a folded skew vector above should not be necessary,
+    #   since the resulting (perhaps) new unit cell does describe the same lattice.
+    #   Therefore the (rotated) atoms are still at correct lattice sites, only some of them
+    #   probably are outside of that new unit cell, which LAMMPS can hopefully deal with.
+    f.write("Atoms \n\n")
+    for i, r in enumerate(atoms.get_positions()):
+        s = species.index(symbols[i]) + 1
+        [x, y, z] = np.dot(r, rotation)
+#        print r, [x, y, z]
+        f.write("%6d %3d %22.15f %22.15f %22.15f \n" % (i+1, s, x, y, z))
+
+    f.close()
+
+
+if __name__ == "__main__":
+
+    pair_style = "eam"
+    Pd_eam_file = "Pd_u3.eam"
+    pair_coeff = [ "* * " + Pd_eam_file ]
+    parameters = { "pair_style" : pair_style, "pair_coeff" : pair_coeff }
+    files = [ Pd_eam_file ]
+    calc = LAMMPS(parameters=parameters, files=files)
+
+    from ase import Atoms
+    a0 = 3.93
+    b0 = a0 / 2.0
+    if True:
+        bulk = Atoms(['Pd']*4,
+                     positions=[(0,0,0),(b0,b0,0),(b0,0,b0),(0,b0,b0)],
+                     cell=[a0]*3,
+                     pbc=True)
+        # test get_forces
+        print "forces for a = %f" % a0
+        print calc.get_forces(bulk)
+        # single points for various lattice constants
+        bulk.set_calculator(calc)
+        for n in range(-5,5,1):
+            a = a0 * (1 + n/100.0)
+            bulk.set_cell([a]*3)
+            print "a : %f , total energy : %f" % (a, bulk.get_potential_energy())
+
+#    calc.clean(force=True)
+    calc.clean()
+
+
+
