|
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>
|