Révision 4 ase/calculators/lammps.py
| lammps.py (revision 4) | ||
|---|---|---|
| 1 |
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> |
|
| 2 |
<html xmlns="http://www.w3.org/1999/xhtml"> |
|
| 3 |
|
|
| 4 |
|
|
| 1 |
#!/usr//bin/env python |
|
| 5 | 2 |
|
| 3 |
# lammps.py (2011/02/22) |
|
| 4 |
# An ASE calculator for the LAMMPS classical MD code available from |
|
| 5 |
# http://lammps.sandia.gov/ |
|
| 6 |
# The environment variable LAMMPS_COMMAND must be defined to point to the LAMMPS binary. |
|
| 7 |
# |
|
| 8 |
# Copyright (C) 2009 - 2011 Joerg Meyer, joerg.meyer@ch.tum.de |
|
| 9 |
# |
|
| 10 |
# This library is free software; you can redistribute it and/or |
|
| 11 |
# modify it under the terms of the GNU Lesser General Public |
|
| 12 |
# License as published by the Free Software Foundation; either |
|
| 13 |
# version 2.1 of the License, or (at your option) any later version. |
|
| 14 |
# |
|
| 15 |
# This library is distributed in the hope that it will be useful, |
|
| 16 |
# but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
| 17 |
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
| 18 |
# Lesser General Public License for more details. |
|
| 19 |
# |
|
| 20 |
# You should have received a copy of the GNU Lesser General Public |
|
| 21 |
# License along with this file; if not, write to the Free Software |
|
| 22 |
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
|
| 23 |
# or see <http://www.gnu.org/licenses/>. |
|
| 6 | 24 |
|
| 7 |
<head> |
|
| 8 |
<title> |
|
| 9 |
lammps.py in trunk/ase/calculators |
|
| 10 |
– ase |
|
| 11 |
</title> |
|
| 12 |
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8" /> |
|
| 13 |
<link rel="search" href="/projects/ase/search" /> |
|
| 14 |
<link rel="help" href="/projects/ase/wiki/TracGuide" /> |
|
| 15 |
<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" /> |
|
| 16 |
<link rel="start" href="/projects/ase/wiki" /> |
|
| 17 |
<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" /> |
|
| 18 |
<link rel="prev" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?rev=2022" title="Revision 2022" /> |
|
| 19 |
<link rel="shortcut icon" href="/projects/ase/chrome/common/trac.ico" type="image/x-icon" /> |
|
| 20 |
<link rel="icon" href="/projects/ase/chrome/common/trac.ico" type="image/x-icon" /> |
|
| 21 |
<link type="application/opensearchdescription+xml" rel="search" href="/projects/ase/search/opensearch" title="Search ase" /> |
|
| 22 |
<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> |
|
| 23 |
<!--[if lt IE 7]> |
|
| 24 |
<script type="text/javascript" src="/projects/ase/chrome/common/js/ie_pre7_hacks.js"></script> |
|
| 25 |
<![endif]--> |
|
| 26 |
<script type="text/javascript" src="/projects/ase/chrome/common/js/folding.js"></script><script type="text/javascript"> |
|
| 27 |
jQuery(document).ready(function($) {
|
|
| 28 |
$(".trac-toggledeleted").show().click(function() {
|
|
| 29 |
$(this).siblings().find(".trac-deleted").toggle();
|
|
| 30 |
return false; |
|
| 31 |
}).click(); |
|
| 32 |
$("#jumploc input").hide();
|
|
| 33 |
$("#jumploc select").change(function () {
|
|
| 34 |
this.parentNode.parentNode.submit(); |
|
| 35 |
}); |
|
| 36 |
$('#preview table.code').enableCollapsibleColumns($('#preview table.code thead th.content'));
|
|
| 37 |
}); |
|
| 38 |
</script> |
|
| 39 |
</head> |
|
| 40 |
<body> |
|
| 41 |
<div id="banner"> |
|
| 42 |
<div id="header"> |
|
| 43 |
<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> |
|
| 44 |
</div> |
|
| 45 |
<form id="search" action="/projects/ase/search" method="get"> |
|
| 46 |
<div> |
|
| 47 |
<label for="proj-search">Search:</label> |
|
| 48 |
<input type="text" id="proj-search" name="q" size="18" value="" /> |
|
| 49 |
<input type="submit" value="Search" /> |
|
| 50 |
</div> |
|
| 51 |
</form> |
|
| 52 |
<div id="metanav" class="nav"> |
|
| 53 |
<ul> |
|
| 54 |
<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> |
|
| 55 |
</ul> |
|
| 56 |
</div> |
|
| 57 |
</div> |
|
| 58 |
<div id="mainnav" class="nav"> |
|
| 59 |
<ul> |
|
| 60 |
<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> |
|
| 61 |
</ul> |
|
| 62 |
</div> |
|
| 63 |
<div id="main"> |
|
| 64 |
<div id="ctxtnav" class="nav"> |
|
| 65 |
<h2>Context Navigation</h2> |
|
| 66 |
<ul> |
|
| 67 |
<li class="first"><span>← <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 →</span></li><li><a href="/projects/ase/browser/trunk/ase/calculators/lammps.py?annotate=blame&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> |
|
| 68 |
</ul> |
|
| 69 |
<hr /> |
|
| 70 |
</div> |
|
| 71 |
<div id="content" class="browser"> |
|
| 72 |
<h1> |
|
| 73 |
<a class="pathentry first" href="/projects/ase/browser?order=name" title="Go to repository root">source:</a> |
|
| 74 |
<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> |
|
| 75 |
<span class="pathentry sep">@</span> |
|
| 76 |
<a class="pathentry" href="/projects/ase/changeset/2023" title="View changeset 2023">2023</a> |
|
| 77 |
<br style="clear: both" /> |
|
| 78 |
</h1> |
|
| 79 |
<div id="jumprev"> |
|
| 80 |
<form action="" method="get"> |
|
| 81 |
<div> |
|
| 82 |
<label for="rev"> |
|
| 83 |
View revision:</label> |
|
| 84 |
<input type="text" id="rev" name="rev" size="6" /> |
|
| 85 |
</div> |
|
| 86 |
</form> |
|
| 87 |
</div> |
|
| 88 |
<div id="jumploc"> |
|
| 89 |
<form action="" method="get"> |
|
| 90 |
<div class="buttons"> |
|
| 91 |
<label for="preselected">Visit:</label> |
|
| 92 |
<select id="preselected" name="preselected"> |
|
| 93 |
<option selected="selected"></option> |
|
| 94 |
<optgroup label="branches"> |
|
| 95 |
<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> |
|
| 96 |
</optgroup><optgroup label="tags"> |
|
| 97 |
<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> |
|
| 98 |
</optgroup> |
|
| 99 |
</select> |
|
| 100 |
<input type="submit" value="Go!" title="Jump to the chosen preselected path" /> |
|
| 101 |
</div> |
|
| 102 |
</form> |
|
| 103 |
</div> |
|
| 104 |
<table id="info" summary="Revision info"> |
|
| 105 |
<tr> |
|
| 106 |
<th scope="col">Revision <a href="/projects/ase/changeset/2023">2023</a>, |
|
| 107 |
<span title="23511 bytes">23.0 KB</span> |
|
| 108 |
checked in by jensj, <a class="timeline" href="/projects/ase/timeline?from=2011-02-23T08%3A39%3A36%2B01%3A00&precision=second" title="2011-02-23T08:39:36+01:00 in Timeline">4 weeks</a> ago |
|
| 109 |
(<a href="/projects/ase/changeset/2023/trunk/ase/calculators/lammps.py">diff</a>)</th> |
|
| 110 |
</tr> |
|
| 111 |
<tr> |
|
| 112 |
<td class="message searchable"> |
|
| 113 |
<p> |
|
| 114 |
Python 2.4 compatibility.<br /> |
|
| 115 |
</p> |
|
| 116 |
</td> |
|
| 117 |
</tr> |
|
| 118 |
</table> |
|
| 119 |
<div id="preview" class="searchable"> |
|
| 25 |
import os |
|
| 26 |
import shutil |
|
| 27 |
import time |
|
| 28 |
import math |
|
| 29 |
import decimal as dec |
|
| 30 |
import numpy as np |
|
| 31 |
from ase import Atoms |
|
| 32 |
from ase.parallel import paropen |
|
| 33 |
|
|
| 34 |
|
|
| 35 |
class LAMMPS: |
|
| 36 |
|
|
| 37 |
def __init__(self, label="lammps", dir="LAMMPS", parameters={}, files=[], always_triclinic=False):
|
|
| 38 |
self.label = label |
|
| 39 |
self.dir = dir |
|
| 40 |
self.parameters = parameters |
|
| 41 |
self.files = files |
|
| 42 |
self.always_triclinic = always_triclinic |
|
| 43 |
self.calls = 0 |
|
| 44 |
self.potential_energy = None |
|
| 45 |
self.forces = None |
|
| 46 |
|
|
| 47 |
if os.path.isdir(self.dir): |
|
| 48 |
if (os.listdir(self.dir) == []): |
|
| 49 |
os.rmdir(self.dir) |
|
| 50 |
else: |
|
| 51 |
os.rename(self.dir, self.dir + ".bak-" + time.strftime("%Y%m%d-%H%M%S"))
|
|
| 52 |
os.mkdir(self.dir, 0755) |
|
| 120 | 53 |
|
| 121 |
<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 <http://www.gnu.org/licenses/>.</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 (>=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) >= 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 >=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) >= 3):</td></tr><tr><th id="L281"><a href="#L281">281</a></th><td> # for >=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) >= 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) <- 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) > 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 -> 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 -> 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> # -> 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 -> 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 -> t_folded=+period/2</td></tr><tr><th id="L495"><a href="#L495">495</a></th><td> # t=+period/2 -> 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> # -> 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>
|
|
| 54 |
for f in files: |
|
| 55 |
shutil.copy(f, os.path.join(self.dir, f)) |
|
| 122 | 56 |
|
| 123 |
</div> |
|
| 124 |
<div id="help"><strong>Note:</strong> See <a href="/projects/ase/wiki/TracBrowser">TracBrowser</a> |
|
| 125 |
for help on using the repository browser.</div> |
|
| 126 |
<div id="anydiff"> |
|
| 127 |
<form action="/projects/ase/diff" method="get"> |
|
| 128 |
<div class="buttons"> |
|
| 129 |
<input type="hidden" name="new_path" value="/trunk/ase/calculators/lammps.py" /> |
|
| 130 |
<input type="hidden" name="old_path" value="/trunk/ase/calculators/lammps.py" /> |
|
| 131 |
<input type="hidden" name="new_rev" /> |
|
| 132 |
<input type="hidden" name="old_rev" /> |
|
| 133 |
<input type="submit" value="View changes..." title="Select paths and revs for Diff" /> |
|
| 134 |
</div> |
|
| 135 |
</form> |
|
| 136 |
</div> |
|
| 137 |
</div> |
|
| 138 |
<div id="altlinks"> |
|
| 139 |
<h3>Download in other formats:</h3> |
|
| 140 |
<ul> |
|
| 141 |
<li class="first"> |
|
| 142 |
<a rel="nofollow" href="/projects/ase/browser/trunk/ase/calculators/lammps.py?format=txt">Plain Text</a> |
|
| 143 |
</li><li class="last"> |
|
| 144 |
<a rel="nofollow" href="/projects/ase/export/2073/trunk/ase/calculators/lammps.py">Original Format</a> |
|
| 145 |
</li> |
|
| 146 |
</ul> |
|
| 147 |
</div> |
|
| 148 |
</div> |
|
| 149 |
<div id="footer" lang="en" xml:lang="en"><hr /> |
|
| 150 |
<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> |
|
| 151 |
<p class="left">Powered by <a href="/projects/ase/about"><strong>Trac 0.12</strong></a><br /> |
|
| 152 |
By <a href="http://www.edgewall.org/">Edgewall Software</a>.</p> |
|
| 153 |
<p class="right">Visit the Trac open source project at<br /><a href="http://trac.edgewall.org/">http://trac.edgewall.org/</a></p> |
|
| 154 |
</div> |
|
| 155 |
</body> |
|
| 156 |
</html> |
|
| 57 |
def clean(self, force=False): |
|
| 58 |
if os.path.isdir(self.dir): |
|
| 59 |
files_in_dir = os.listdir(self.dir) |
|
| 60 |
files_copied = self.files |
|
| 61 |
if (len(files_in_dir) == 0): |
|
| 62 |
os.rmdir(self.dir) |
|
| 63 |
elif (set(files_in_dir).issubset(set(files_copied)) or force): |
|
| 64 |
shutil.rmtree(self.dir) |
|
| 65 |
|
|
| 66 |
def get_potential_energy(self, atoms): |
|
| 67 |
self.update(atoms) |
|
| 68 |
return self.potential_energy |
|
| 69 |
|
|
| 70 |
def get_forces(self, atoms): |
|
| 71 |
self.update(atoms) |
|
| 72 |
return self.forces.copy() |
|
| 73 |
|
|
| 74 |
def get_stress(self, atoms): |
|
| 75 |
raise NotImplementedError |
|
| 76 |
|
|
| 77 |
def update(self, atoms): |
|
| 78 |
# TODO: check if (re-)calculation is necessary |
|
| 79 |
self.calculate(atoms) |
|
| 80 |
|
|
| 81 |
def calculate(self, atoms): |
|
| 82 |
self.atoms = atoms.copy() |
|
| 83 |
self.run() |
|
| 84 |
|
|
| 85 |
def run(self): |
|
| 86 |
"""Method which explicitely runs LAMMPS.""" |
|
| 87 |
|
|
| 88 |
self.calls += 1 |
|
| 89 |
|
|
| 90 |
# set LAMMPS command from environment variable |
|
| 91 |
if os.environ.has_key('LAMMPS_COMMAND'):
|
|
| 92 |
lammps_command = os.path.abspath(os.environ['LAMMPS_COMMAND']) |
|
| 93 |
else: |
|
| 94 |
self.clean() |
|
| 95 |
raise RuntimeError('Please set LAMMPS_COMMAND environment variable')
|
|
| 96 |
if os.environ.has_key('LAMMPS_OPTIONS'):
|
|
| 97 |
lammps_options = os.environ['LAMMPS_OPTIONS'] |
|
| 98 |
else: |
|
| 99 |
lammps_options = "-echo log -screen none" |
|
| 100 |
|
|
| 101 |
# change into subdirectory for LAMMPS calculations |
|
| 102 |
cwd = os.getcwd() |
|
| 103 |
os.chdir(self.dir) |
|
| 104 |
|
|
| 105 |
# setup file names for LAMMPS calculation |
|
| 106 |
label = "%s-%06d" % (self.label, self.calls) |
|
| 107 |
lammps_data = "data." + label |
|
| 108 |
lammps_in = "in." + label |
|
| 109 |
lammps_trj = label + ".lammpstrj" |
|
| 110 |
lammps_log = label + ".log" |
|
| 111 |
|
|
| 112 |
# write LAMMPS input files |
|
| 113 |
self.write_lammps_data(lammps_data=lammps_data) |
|
| 114 |
self.write_lammps_in(lammps_in=lammps_in, lammps_data=lammps_data, lammps_trj=lammps_trj, parameters=self.parameters) |
|
| 115 |
|
|
| 116 |
# run LAMMPS |
|
| 117 |
# TODO: check for successful completion (based on output files?!) of LAMMPS call... |
|
| 118 |
exitcode = os.system("%s %s -in %s -log %s" % (lammps_command, lammps_options, lammps_in, lammps_log))
|
|
| 119 |
if exitcode != 0: |
|
| 120 |
cwd = os.getcwd() |
|
| 121 |
raise RuntimeError("LAMMPS exited in %s with exit code: %d. " % (cwd,exitcode))
|
|
| 122 |
|
|
| 123 |
# read LAMMPS output files |
|
| 124 |
self.read_lammps_log(lammps_log=lammps_log) |
|
| 125 |
self.read_lammps_trj(lammps_trj=lammps_trj) |
|
| 126 |
|
|
| 127 |
# change back to previous working directory |
|
| 128 |
os.chdir(cwd) |
|
| 129 |
|
|
| 130 |
def write_lammps_data(self, lammps_data=None): |
|
| 131 |
"""Method which writes a LAMMPS data file with atomic structure.""" |
|
| 132 |
if (lammps_data == None): |
|
| 133 |
lammps_data = "data." + self.label |
|
| 134 |
write_lammps(lammps_data, self.atoms, force_skew=self.always_triclinic) |
|
| 135 |
|
|
| 136 |
def write_lammps_in(self, lammps_in=None, lammps_data=None, lammps_trj=None, parameters={}):
|
|
| 137 |
"""Method which writes a LAMMPS in file with run parameters and settings.""" |
|
| 138 |
if (lammps_in == None): |
|
| 139 |
lammps_in = "in." + self.label |
|
| 140 |
if (lammps_data == None): |
|
| 141 |
lammps_data = "data." + self.label |
|
| 142 |
if (lammps_trj == None): |
|
| 143 |
lammps_trj = self.label + ".lammpstrj" |
|
| 144 |
|
|
| 145 |
f = paropen(lammps_in, 'w') |
|
| 146 |
f.write("# " + f.name + " (written by ASE) \n")
|
|
| 147 |
f.write("\n")
|
|
| 148 |
|
|
| 149 |
f.write("### variables \n")
|
|
| 150 |
f.write("variable data_file index \"%s\" \n" % lammps_data)
|
|
| 151 |
f.write("variable dump_file index \"%s\" \n" % lammps_trj)
|
|
| 152 |
f.write("\n\n")
|
|
| 153 |
|
|
| 154 |
pbc = self.atoms.get_pbc() |
|
| 155 |
f.write("### simulation box \n")
|
|
| 156 |
f.write("units metal \n")
|
|
| 157 |
if ("boundary" in parameters):
|
|
| 158 |
f.write("boundary %s \n" % parameters["boundary"])
|
|
| 159 |
else: |
|
| 160 |
f.write("boundary %c %c %c \n" % tuple('sp'[x] for x in pbc))
|
|
| 161 |
f.write("atom_modify sort 0 0.0 \n")
|
|
| 162 |
for key in ("neighbor" ,"newton"):
|
|
| 163 |
if key in parameters: |
|
| 164 |
f.write("%s %s \n" % (key, parameters[key]))
|
|
| 165 |
f.write("\n")
|
|
| 166 |
f.write("read_data ${data_file} \n")
|
|
| 167 |
f.write("\n\n")
|
|
| 168 |
|
|
| 169 |
f.write("### interactions \n")
|
|
| 170 |
if ( ("pair_style" in parameters) and ("pair_coeff" in parameters) ):
|
|
| 171 |
pair_style = parameters["pair_style"] |
|
| 172 |
f.write("pair_style %s \n" % pair_style)
|
|
| 173 |
for pair_coeff in parameters["pair_coeff"]: |
|
| 174 |
f.write("pair_coeff %s \n" % pair_coeff)
|
|
| 175 |
if "mass" in parameters: |
|
| 176 |
for mass in parameters["mass"]: |
|
| 177 |
f.write("mass %s \n" % mass)
|
|
| 178 |
else: |
|
| 179 |
# simple default parameters that should always make the LAMMPS calculation run |
|
| 180 |
# pair_style lj/cut 2.5 |
|
| 181 |
# pair_coeff * * 1 1 |
|
| 182 |
# mass * 1.0 else: |
|
| 183 |
f.write("pair_style lj/cut 2.5 \n")
|
|
| 184 |
f.write("pair_coeff * * 1 1 \n")
|
|
| 185 |
f.write("mass * 1.0 \n")
|
|
| 186 |
f.write("\n")
|
|
| 187 |
|
|
| 188 |
f.write("### run \n")
|
|
| 189 |
f.write("fix fix_nve all nve \n")
|
|
| 190 |
f.write("\n")
|
|
| 191 |
f.write("dump dump_all all custom 1 ${dump_file} id type x y z vx vy vz fx fy fz \n")
|
|
| 192 |
f.write("\n")
|
|
| 193 |
f.write("thermo_style custom step temp ke pe etotal cpu \n")
|
|
| 194 |
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")
|
|
| 195 |
f.write("thermo 1 \n")
|
|
| 196 |
f.write("\n")
|
|
| 197 |
if ("minimize" in parameters):
|
|
| 198 |
f.write("minimize %s \n" % parameters["minimize"])
|
|
| 199 |
if ("run" in parameters):
|
|
| 200 |
f.write("run %s \n" % parameters["run"])
|
|
| 201 |
if not ( ("minimize" in parameters) or ("run" in parameters) ):
|
|
| 202 |
f.write("run 0 \n")
|
|
| 203 |
|
|
| 204 |
f.close() |
|
| 205 |
|
|
| 206 |
def read_lammps_log(self, lammps_log=None): |
|
| 207 |
"""Method which reads a LAMMPS output log file.""" |
|
| 208 |
if (lammps_log == None): |
|
| 209 |
lammps_log = self.label + ".log" |
|
| 210 |
|
|
| 211 |
epot = 0.0 |
|
| 212 |
|
|
| 213 |
f = paropen(lammps_log, 'r') |
|
| 214 |
while True: |
|
| 215 |
line = f.readline() |
|
| 216 |
if not line: |
|
| 217 |
break |
|
| 218 |
# get potential energy of first step (if more are done) |
|
| 219 |
if "PotEng" in line: |
|
| 220 |
i = line.split().index("PotEng")
|
|
| 221 |
line = f.readline() |
|
| 222 |
epot = float(line.split()[i]) |
|
| 223 |
f.close() |
|
| 224 |
|
|
| 225 |
# print epot |
|
| 226 |
|
|
| 227 |
self.potential_energy = epot |
|
| 228 |
|
|
| 229 |
def read_lammps_trj(self, lammps_trj=None, set_atoms=False): |
|
| 230 |
"""Method which reads a LAMMPS dump file.""" |
|
| 231 |
if (lammps_trj == None): |
|
| 232 |
lammps_trj = self.label + ".lammpstrj" |
|
| 233 |
|
|
| 234 |
f = paropen(lammps_trj, 'r') |
|
| 235 |
while True: |
|
| 236 |
line = f.readline() |
|
| 237 |
|
|
| 238 |
if not line: |
|
| 239 |
break |
|
| 240 |
|
|
| 241 |
#TODO: extend to proper dealing with multiple steps in one trajectory file |
|
| 242 |
if "ITEM: TIMESTEP" in line: |
|
| 243 |
n_atoms = 0 |
|
| 244 |
lo = [] ; hi = [] ; tilt = [] |
|
| 245 |
id = [] ; type = [] |
|
| 246 |
positions = [] ; velocities = [] ; forces = [] |
|
| 247 |
|
|
| 248 |
if "ITEM: NUMBER OF ATOMS" in line: |
|
| 249 |
line = f.readline() |
|
| 250 |
n_atoms = int(line.split()[0]) |
|
| 251 |
|
|
| 252 |
if "ITEM: BOX BOUNDS" in line: |
|
| 253 |
# save labels behind "ITEM: BOX BOUNDS" in triclinic case (>=lammps-7Jul09) |
|
| 254 |
tilt_items = line.split()[3:] |
|
| 255 |
for i in range(3): |
|
| 256 |
line = f.readline() |
|
| 257 |
fields = line.split() |
|
| 258 |
lo.append(float(fields[0])) |
|
| 259 |
hi.append(float(fields[1])) |
|
| 260 |
if (len(fields) >= 3): |
|
| 261 |
tilt.append(float(fields[2])) |
|
| 262 |
|
|
| 263 |
if "ITEM: ATOMS" in line: |
|
| 264 |
# (reliably) identify values by labels behind "ITEM: ATOMS" - requires >=lammps-7Jul09 |
|
| 265 |
# create corresponding index dictionary before iterating over atoms to (hopefully) speed up lookups... |
|
| 266 |
atom_attributes = {}
|
|
| 267 |
for (i, x) in enumerate(line.split()[2:]): |
|
| 268 |
atom_attributes[x] = i |
|
| 269 |
for n in range(n_atoms): |
|
| 270 |
line = f.readline() |
|
| 271 |
fields = line.split() |
|
| 272 |
id.append( int(fields[atom_attributes['id']]) ) |
|
| 273 |
type.append( int(fields[atom_attributes['type']]) ) |
|
| 274 |
positions.append( [ float(fields[atom_attributes[x]]) for x in ['x', 'y', 'z'] ] ) |
|
| 275 |
velocities.append( [ float(fields[atom_attributes[x]]) for x in ['vx', 'vy', 'vz'] ] ) |
|
| 276 |
forces.append( [ float(fields[atom_attributes[x]]) for x in ['fx', 'fy', 'fz'] ] ) |
|
| 277 |
f.close() |
|
| 278 |
|
|
| 279 |
# determine cell tilt (triclinic case!) |
|
| 280 |
if (len(tilt) >= 3): |
|
| 281 |
# for >=lammps-7Jul09 use labels behind "ITEM: BOX BOUNDS" to assign tilt (vector) elements ... |
|
| 282 |
if (len(tilt_items) >= 3): |
|
| 283 |
xy = tilt[tilt_items.index('xy')]
|
|
| 284 |
xz = tilt[tilt_items.index('xz')]
|
|
| 285 |
yz = tilt[tilt_items.index('yz')]
|
|
| 286 |
# ... otherwise assume default order in 3rd column (if the latter was present) |
|
| 287 |
else: |
|
| 288 |
xy = tilt[0] |
|
| 289 |
xz = tilt[1] |
|
| 290 |
yz = tilt[2] |
|
| 291 |
else: |
|
| 292 |
xy = xz = yz = 0 |
|
| 293 |
xhilo = (hi[0] - lo[0]) - xy - xz |
|
| 294 |
yhilo = (hi[1] - lo[1]) - yz |
|
| 295 |
zhilo = (hi[2] - lo[2]) |
|
| 296 |
|
|
| 297 |
# The simulation box bounds are included in each snapshot and if the box is triclinic (non-orthogonal), |
|
| 298 |
# then the tilt factors are also printed; see the region prism command for a description of tilt factors. |
|
| 299 |
# For triclinic boxes the box bounds themselves (first 2 quantities on each line) are a true "bounding box" |
|
| 300 |
# around the simulation domain, which means they include the effect of any tilt. |
|
| 301 |
# [ http://lammps.sandia.gov/doc/dump.html , lammps-7Jul09 ] |
|
| 302 |
# |
|
| 303 |
# This *should* extract the lattice vectors that LAMMPS uses from the true "bounding box" printed in the dump file |
|
| 304 |
# It might fail in some cases (negative tilts?!) due to the MIN / MAX construction of these box corners: |
|
| 305 |
# |
|
| 306 |
# void Domain::set_global_box() |
|
| 307 |
# [...] |
|
| 308 |
# if (triclinic) {
|
|
| 309 |
# [...] |
|
| 310 |
# boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy); |
|
| 311 |
# boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz); |
|
| 312 |
# boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz); |
|
| 313 |
# boxlo_bound[2] = boxlo[2]; |
|
| 314 |
# |
|
| 315 |
# boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy); |
|
| 316 |
# boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz); |
|
| 317 |
# boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz); |
|
| 318 |
# boxhi_bound[2] = boxhi[2]; |
|
| 319 |
# } |
|
| 320 |
# [ lammps-7Jul09/src/domain.cpp ] |
|
| 321 |
# |
|
| 322 |
cell = [[xhilo,0,0],[xy,yhilo,0],[xz,yz,zhilo]] |
|
| 323 |
|
|
| 324 |
# print n_atoms |
|
| 325 |
# print lo |
|
| 326 |
# print hi |
|
| 327 |
# print id |
|
| 328 |
# print type |
|
| 329 |
# print positions |
|
| 330 |
# print velocities |
|
| 331 |
# print forces |
|
| 332 |
|
|
| 333 |
# assume that LAMMPS does not reorder atoms internally |
|
| 334 |
# print np.array(positions) |
|
| 335 |
# print self.atoms.get_positions() |
|
| 336 |
# print np.array(positions) - self.atoms.get_positions() |
|
| 337 |
|
|
| 338 |
cell_atoms = np.array(cell) |
|
| 339 |
type_atoms = np.array(type) |
|
| 340 |
# velocities_atoms = np.array(velocities) |
|
| 341 |
positions_atoms = np.array(positions) |
|
| 342 |
forces_atoms = np.array(forces) |
|
| 343 |
|
|
| 344 |
if self.atoms: |
|
| 345 |
cell_atoms = self.atoms.get_cell() |
|
| 346 |
|
|
| 347 |
rotation_lammps2ase = np.dot(np.linalg.inv(np.array(cell)), cell_atoms) |
|
| 348 |
# print "rotation_lammps2ase:" |
|
| 349 |
# print rotation_lammps2ase |
|
| 350 |
# print np.transpose(rotation_lammps2ase) |
|
| 351 |
# print np.linalg.det(rotation_lammps2ase) # check for orthogonality of matrix 'rotation' |
|
| 352 |
# print np.dot(rotation_lammps2ase, np.transpose(rotation_lammps2ase)) # check for orthogonality of matrix 'rotation' |
|
| 353 |
|
|
| 354 |
type_atoms = self.atoms.get_atomic_numbers() |
|
| 355 |
positions_atoms = np.array( [np.dot(np.array(r), rotation_lammps2ase) for r in positions] ) |
|
| 356 |
velocities_atoms = np.array( [np.dot(np.array(v), rotation_lammps2ase) for v in velocities] ) |
|
| 357 |
forces_atoms = np.array( [np.dot(np.array(f), rotation_lammps2ase) for f in forces] ) |
|
| 358 |
|
|
| 359 |
if (set_atoms): |
|
| 360 |
# assume periodic boundary conditions here (like also below in write_lammps) <- TODO:?! |
|
| 361 |
self.atoms = Atoms(type_atoms, positions=positions_atoms, cell=cell_atoms, pbc=True) |
|
| 362 |
|
|
| 363 |
self.forces = forces_atoms |
|
| 364 |
|
|
| 365 |
|
|
| 366 |
# could perhaps go into io.lammps in the future... |
|
| 367 |
#from ase.parallel import paropen |
|
| 368 |
def write_lammps(fileobj, atoms, force_skew=False): |
|
| 369 |
"""Method which writes atomic structure data to a LAMMPS data file.""" |
|
| 370 |
if isinstance(fileobj, file): |
|
| 371 |
f = fileobj |
|
| 372 |
elif isinstance(fileobj, str): |
|
| 373 |
f = paropen(fileobj, 'w') |
|
| 374 |
else : |
|
| 375 |
raise TypeError('fileobj is not of type file or str!')
|
|
| 376 |
|
|
| 377 |
if isinstance(atoms, list): |
|
| 378 |
if len(atoms) > 1: |
|
| 379 |
raise ValueError('Can only write one configuration to a lammps data file!')
|
|
| 380 |
atoms = atoms[0] |
|
| 381 |
|
|
| 382 |
f.write(f.name + " (written by ASE) \n\n") |
|
| 383 |
|
|
| 384 |
symbols = atoms.get_chemical_symbols() |
|
| 385 |
n_atoms = len(symbols) |
|
| 386 |
f.write("%d \t atoms \n" % n_atoms)
|
|
| 387 |
|
|
| 388 |
# Uniqify 'symbols' list and sort to a new list 'species'. |
|
| 389 |
# Sorting is important in case of multi-component systems: |
|
| 390 |
# This way it is assured that LAMMPS atom types are always |
|
| 391 |
# assigned predictively according to the alphabetic order |
|
| 392 |
# of the species present in the current system. |
|
| 393 |
# (Hence e.g. the mapping in interaction statements for alloy |
|
| 394 |
# potentials depending on the order of the elements in the |
|
| 395 |
# external potential can also be safely adjusted accordingly.) |
|
| 396 |
species = sorted(list(set(symbols))) |
|
| 397 |
n_atom_types = len(species) |
|
| 398 |
f.write("%d \t atom types \n" % n_atom_types)
|
|
| 399 |
|
|
| 400 |
|
|
| 401 |
### cell -> bounding box for LAMMPS |
|
| 402 |
|
|
| 403 |
## A) simple version |
|
| 404 |
# - For orthogonal cells which are properly aligned with the reference coordinate system only. |
|
| 405 |
# cell = atoms.get_cell() |
|
| 406 |
# lo = [0.0, 0.0, 0.0] |
|
| 407 |
# hi = sum(cell[0:3]) |
|
| 408 |
# f.write("%15.8f %15.8f \t xlo xhi \n" % (lo[0], hi[0]))
|
|
| 409 |
# f.write("%15.8f %15.8f \t ylo yhi \n" % (lo[1], hi[1]))
|
|
| 410 |
# f.write("%15.8f %15.8f \t zlo zhi \n" % (lo[2], hi[2]))
|
|
| 411 |
# f.write("\n\n")
|
|
| 412 |
|
|
| 413 |
## B) convert to (perhaps) triclinic cells in LAMMPS |
|
| 414 |
|
|
| 415 |
# B.1) calculate lengths of and angles between cell vectors |
|
| 416 |
# (conventional crystallographic nomenclature) |
|
| 417 |
cell = atoms.get_cell() |
|
| 418 |
a = np.linalg.norm(cell[0]) |
|
| 419 |
b = np.linalg.norm(cell[1]) |
|
| 420 |
c = np.linalg.norm(cell[2]) |
|
| 421 |
alpha = np.arccos( np.vdot(cell[1], cell[2]) / (b * c) ) |
|
| 422 |
beta = np.arccos( np.vdot(cell[0], cell[2]) / (a * c) ) |
|
| 423 |
gamma = np.arccos( np.vdot(cell[0], cell[1]) / (a * b) ) |
|
| 424 |
# print [a, b, c] |
|
| 425 |
# print map(np.degrees, [alpha, beta, gamma]) |
|
| 426 |
|
|
| 427 |
# B.2) construct bounding box edges and skew vector of |
|
| 428 |
# corresponding (perhaps rotated!) LAMMPS cell |
|
| 429 |
# http://lammps.sandia.gov/doc/read_data.html : |
|
| 430 |
# a_LAMMPS = (xhi-xlo,0,0); b_LAMMPS = (xy,yhi-ylo,0); c_LAMMPS = (xz,yz,zhi-zlo) |
|
| 431 |
xlo = ylo = zlo = 0.0 |
|
| 432 |
# this choice of origin simplifies things: |
|
| 433 |
# a_LAMMPS = (xhi,0,0); b_LAMMPS = (xy,yhi,0); c_LAMMPS = (xz,yz,zhi) |
|
| 434 |
xhi = a # a_LAMMPS |
|
| 435 |
xy = np.cos(gamma) * b # b_LAMMPS |
|
| 436 |
yhi = np.sin(gamma) * b |
|
| 437 |
xz = np.cos(beta) * c # c_LAMMPS |
|
| 438 |
yz = ( b * c * np.cos(alpha) - xy * xz ) / yhi |
|
| 439 |
zhi = np.sqrt( c**2 - xz**2 - yz**2 ) |
|
| 440 |
|
|
| 441 |
# B.3) obtain rotation of original cell with respect to LAMMPS cell |
|
| 442 |
# to properly tranform the atoms contained within (see below!) |
|
| 443 |
# IMPORTANT: use same convention as in cell of ASE atoms object, i.e. |
|
| 444 |
# cell vectors are ROW VECTORS in numpy array |
|
| 445 |
cell_lammps = np.array([[xhi-xlo,0,0],[xy,yhi-ylo,0],[xz,yz,zhi-zlo]]) |
|
| 446 |
# a_lammps = np.linalg.norm(cell_lammps[0]) |
|
| 447 |
# b_lammps = np.linalg.norm(cell_lammps[1]) |
|
| 448 |
# c_lammps = np.linalg.norm(cell_lammps[2]) |
|
| 449 |
# alpha_lammps = np.arccos( np.vdot(cell_lammps[1], cell_lammps[2]) / (b_lammps * c_lammps) ) |
|
| 450 |
# beta_lammps = np.arccos( np.vdot(cell_lammps[0], cell_lammps[2]) / (a_lammps * c_lammps) ) |
|
| 451 |
# gamma_lammps = np.arccos( np.vdot(cell_lammps[0], cell_lammps[1]) / (a_lammps * b_lammps) ) |
|
| 452 |
# print [a_lammps, b_lammps, c_lammps] |
|
| 453 |
# print map(np.degrees, [alpha_lammps, beta_lammps, gamma_lammps]) |
|
| 454 |
# IMPORTANT: need vector-(rotation-)matrix product (instead of matrix-vector product) here, |
|
| 455 |
# cell vectors are ROW VECTORS (also see above) |
|
| 456 |
rotation = np.dot(np.linalg.inv(cell), cell_lammps) |
|
| 457 |
# print rotation |
|
| 458 |
# print np.transpose(rotation) |
|
| 459 |
# print np.linalg.det(rotation) # check for orthogonality of matrix 'rotation' |
|
| 460 |
# print np.dot(rotation, np.transpose(rotation)) # check for orthogonality of matrix 'rotation' |
|
| 461 |
# print np.dot(rotation, cell[0]) |
|
| 462 |
# print np.dot(rotation, cell[1]) |
|
| 463 |
# print np.dot(rotation, cell[2]) |
|
| 464 |
|
|
| 465 |
# B.4.1) write bounding box edges |
|
| 466 |
f.write("%15.8f %15.8f \t\t\t xlo xhi \n" % (xlo, xhi))
|
|
| 467 |
f.write("%15.8f %15.8f \t\t\t ylo yhi \n" % (ylo, yhi))
|
|
| 468 |
f.write("%15.8f %15.8f \t\t\t zlo zhi \n" % (zlo, zhi))
|
|
| 469 |
|
|
| 470 |
# B.4.2) sanitize and write skew vector (if necessary) |
|
| 471 |
# The too simple version |
|
| 472 |
# f.write("%20.10f %20.10f %20.10f \t xy xz yz \n" % (xy, xz, yz))
|
|
| 473 |
# can make LAMMPS (easily) reject the definition of a triclinic box because |
|
| 474 |
# a) skew vector components calculated above are outside the respective ranges |
|
| 475 |
# which LAMMPS expects as described here |
|
| 476 |
# http://lammps.sandia.gov/doc/read_data.html |
|
| 477 |
# and coded here |
|
| 478 |
# domain.cpp -> Domain::set_initial_box() |
|
| 479 |
# b) rounding up in the skew vector output can make the result violate those |
|
| 480 |
# ranges for purely numeric reasons (perhaps even more frequent problem than a) |
|
| 481 |
# and definitely more stupid!) |
|
| 482 |
# More elaborate solution addressing the issues described above: |
|
| 483 |
# -> a): |
|
| 484 |
# Fold back skew vector components calculated above to equivalent ones complying with |
|
| 485 |
# what LAMMPS expects (see references above) and describing another ("less skewed")
|
|
| 486 |
# unit cell for the same lattice: |
|
| 487 |
# t -> t_folded in [-period/2, +period/2] |
|
| 488 |
def fold_skew_dec(t, period) : |
|
| 489 |
t_dec = dec.Decimal(repr(t)) |
|
| 490 |
period_dec = dec.Decimal(repr(period)) |
|
| 491 |
# dec.ROUND_HALF_DOWN leaves all t in [-period/2, +period/2] unchanged whereas |
|
| 492 |
# dec.ROUND_HALF_UP does the same for all t in (-period/2, +period/2) but |
|
| 493 |
# swaps the boundary values: |
|
| 494 |
# t=-period/2 -> t_folded=+period/2 |
|
| 495 |
# t=+period/2 -> t_folded=-period/2 |
|
| 496 |
t_folded_dec = t_dec - period_dec * (t_dec/period_dec).to_integral_value(dec.ROUND_HALF_DOWN) |
|
| 497 |
return t_folded_dec |
|
| 498 |
skew_folded_dec = [ fold_skew_dec(xy, xhi-xlo), fold_skew_dec(xz, xhi-xlo), fold_skew_dec(yz, yhi-ylo) ] |
|
| 499 |
# -> b): |
|
| 500 |
# This should be kept in sync with the accuracy used for writing the bounding box edges above |
|
| 501 |
prec_dec = dec.Decimal('1E-8')
|
|
| 502 |
# Make sure to always round down and hence avoid numerical problems with skew vector in LAMMPS. |
|
| 503 |
skew_folded_and_rounded_dec = [ d.quantize(prec_dec, dec.ROUND_DOWN) for d in skew_folded_dec ] |
|
| 504 |
# Do not write zero skew vector. |
|
| 505 |
# (I.e. calculate orthogonal cell with LAMMPS in that case!) |
|
| 506 |
has_skew = any( [not d.is_zero() for d in skew_folded_and_rounded_dec] ) |
|
| 507 |
if force_skew or has_skew: |
|
| 508 |
f.write( "%15s %15s %15s \t xy xz yz \n" % tuple( str(d) for d in skew_folded_and_rounded_dec ) ) |
|
| 509 |
f.write("\n\n")
|
|
| 510 |
|
|
| 511 |
|
|
| 512 |
### atoms |
|
| 513 |
|
|
| 514 |
## A) simple version |
|
| 515 |
# - For orthogonal cells which are properly aligned with the reference coordinate system only. |
|
| 516 |
# f.write("Atoms \n\n")
|
|
| 517 |
# for i, (x, y, z) in enumerate(atoms.get_positions()): |
|
| 518 |
# s = species.index(symbols[i]) + 1 |
|
| 519 |
# f.write("%6d %3d %22.15f %22.15f %22.15f \n" % (i+1, s, x, y, z))
|
|
| 520 |
|
|
| 521 |
## B) convert to (perhaps) triclinic cells in LAMMPS: |
|
| 522 |
# - Apply rotation determined above for this case. |
|
| 523 |
# - Lattice translation in case of a folded skew vector above should not be necessary, |
|
| 524 |
# since the resulting (perhaps) new unit cell does describe the same lattice. |
|
| 525 |
# Therefore the (rotated) atoms are still at correct lattice sites, only some of them |
|
| 526 |
# probably are outside of that new unit cell, which LAMMPS can hopefully deal with. |
|
| 527 |
f.write("Atoms \n\n")
|
|
| 528 |
for i, r in enumerate(atoms.get_positions()): |
|
| 529 |
s = species.index(symbols[i]) + 1 |
|
| 530 |
[x, y, z] = np.dot(r, rotation) |
|
| 531 |
# print r, [x, y, z] |
|
| 532 |
f.write("%6d %3d %22.15f %22.15f %22.15f \n" % (i+1, s, x, y, z))
|
|
| 533 |
|
|
| 534 |
f.close() |
|
| 535 |
|
|
| 536 |
|
|
| 537 |
if __name__ == "__main__": |
|
| 538 |
|
|
| 539 |
pair_style = "eam" |
|
| 540 |
Pd_eam_file = "Pd_u3.eam" |
|
| 541 |
pair_coeff = [ "* * " + Pd_eam_file ] |
|
| 542 |
parameters = { "pair_style" : pair_style, "pair_coeff" : pair_coeff }
|
|
| 543 |
files = [ Pd_eam_file ] |
|
| 544 |
calc = LAMMPS(parameters=parameters, files=files) |
|
| 545 |
|
|
| 546 |
from ase import Atoms |
|
| 547 |
a0 = 3.93 |
|
| 548 |
b0 = a0 / 2.0 |
|
| 549 |
if True: |
|
| 550 |
bulk = Atoms(['Pd']*4, |
|
| 551 |
positions=[(0,0,0),(b0,b0,0),(b0,0,b0),(0,b0,b0)], |
|
| 552 |
cell=[a0]*3, |
|
| 553 |
pbc=True) |
|
| 554 |
# test get_forces |
|
| 555 |
print "forces for a = %f" % a0 |
|
| 556 |
print calc.get_forces(bulk) |
|
| 557 |
# single points for various lattice constants |
|
| 558 |
bulk.set_calculator(calc) |
|
| 559 |
for n in range(-5,5,1): |
|
| 560 |
a = a0 * (1 + n/100.0) |
|
| 561 |
bulk.set_cell([a]*3) |
|
| 562 |
print "a : %f , total energy : %f" % (a, bulk.get_potential_energy()) |
|
| 563 |
|
|
| 564 |
# calc.clean(force=True) |
|
| 565 |
calc.clean() |
|
| 566 |
|
|
| 567 |
|
|
| 568 |
|
|
Formats disponibles : Unified diff