Volume optimization of Ag

Volume optimization is a simple case of an aggregation property. This example shows how to set up the inputs for the data points with the help of a template and how to determine the optimal volume by fitting the corresponding total energies.

We are going to use a simple inputxml template to generate the input files for the total-energies calculations.

Experiment description

In order to vary the scale of the unit cell by +/- 5% we need a description of our experiment:

File: scale.xml

<?xml version="1.0" encoding="UTF-8" ?>
 <experiment>
    <set unicellscaling="0.95" path="unicellscaling_0.95"/>
    <set unicellscaling="0.975" path="unicellscaling_0.975" />
    <set unicellscaling="1" path="unicellscaling_1"/>
    <set unicellscaling="1.025" path="unicellscaling_1.025"/>
    <set unicellscaling="1.05" path="unicellscaling_1.05"/>
 </experiment>

This quasi format is referred to as experiment-list and used in various other examples.

The template

In order to generate the input files we use XSLT templating-tools. The corresponding template:

File:evTemplate.xsl

<?xml version="1.0" encoding="UTF-8" ?>
<xsl:stylesheet version="1.0"
  xmlns:xsl="http://www.w3.org/1999/XSL/Transform">
<xsl:output method="xml" />
<xsl:template match="/">
 
<!-- Authors: chm, jus, sag -->
 
<!-- Loop over all elements named "set" from reference xml-file -->
<xsl:for-each select = "/experiment/set">
<xsl:variable name="path">
<xsl:value-of select="@path"/>
<xsl:text>/imput.xml</xsl:text>
</xsl:variable> 
 
 <!-- Write document at Path $path -->
 <xsl:document href="{$path}" method="xml" indent="yes">
 <xsl:comment>
 This file is generated with XSLTPROC using a template file and a reference file
</xsl:comment>
<!-- ////////////////////////////////////////////////////////////////////////-->
<!-- /// The input file begins here /////////////////////////////////////////-->
<!-- ////////////////////////////////////////////////////////////////////////-->
<input xsi:noNamespaceSchemaLocation="excitinginput.xsd"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsltpath="../../../xml/"
  scratchpath="/tmp/chm/1">
  <title>Ag Energy-Volume</title>
  <structure speciespath="../../../species/">
    <crystal >
    <xsl:attribute name="scale">
    <xsl:value-of select="3.86*@unicellscaling"/>
    </xsl:attribute>
 
      <basevect>1.0 1.0 0.0</basevect>
      <basevect>1.0 0.0 1.0</basevect>
      <basevect>0.0 1.0 1.0</basevect>
    </crystal>
    <species speciesfile="Ag.xml">
      <atom coord="0.0  0.0  0.0" />
    </species>
  </structure>
  <groundstate ngridk="8 8 8" rgkmax="8" />
</input>
 
<!-- ////////////////////////////////////////////////////////////////////////-->
<!-- /// The input file ends here ///////////////////////////////////////////-->
<!-- ////////////////////////////////////////////////////////////////////////-->
 
</xsl:document>
</xsl:for-each>
</xsl:template>
</xsl:stylesheet>

It looks quite complicated, but it consists only of a generic part (which can be used for all cases with a similar experiment description) and the Ag input file, where in the latter the @scale attribute of the <crystal> element is expressed in XSLT:

<crystal>
    <xsl:attribute name="scale">
    <xsl:value-of select="3.86*@unicellscaling"/>
    </xsl:attribute>
    [...]
</crystal>

To apply this template, execute:

xsltproc  evTemplate.xsl  scales.xml

Then xsltproc will create the five directories with the proper input file therein.

Execution

Here we have to run 5 self-consistent calculations. You may start them by hand. For more complex examples one could use templates to create scripts of all kind like: loadleveler-multistep-job-template

Analysis

In order to get the energy-volume curve we have to collect the data, which are the volumes and corresponding total-energy values. They are saved among many other values in info.xml. To collect them into one file we use getev.xsl to get:

<?xml version="1.0"?>
<graph xmlns:math="http://exslt.org/math" xmlns:str="http://exslt.org/strings">
  <point volume="139.469011682385" scale="3.667" totalEnergy="-5312.439270664138"/>
  <point volume="150.772020410909" scale="3.7635" totalEnergy="-5312.4416556378"/>
  <point volume="162.669790561172" scale="3.86" totalEnergy="-5312.440895946467"/>
  <point volume="175.177572426039" scale="3.9565" totalEnergy="-5312.437999079701"/>
  <point volume="188.310616298377" scale="4.053" totalEnergy="-5312.433397077082"/>
</graph>

xsltproc getev.xsl scales.xml >eV.xml

The script EV.py fits the data points, gets the minimum and plots the corresponding graphics:

import numpy as nm
import libxml2
from array import *
from libxml2 import xmlAttr
import matplotlib.pyplot as plt
# read data
eVdoc = libxml2.parseFile("./eV.xml")
ctxt = eVdoc.xpathNewContext()
Volume=nm.array(map(float,map(xmlAttr.getContent,ctxt.xpathEval("//@volume"))))
totalEnergy=nm.array(map(float,map(xmlAttr.getContent,ctxt.xpathEval("//@totalEnergy"))))
# make parabular fit
p=nm.polyfit(Volume,totalEnergy,2)
curve=nm.poly1d(p)
# find root of derivative to get minimum
min=nm.roots(nm.polyder(p))
# x values for plotting polynomial
xa = nm.linspace(Volume[0],Volume[-1],100)
#plot
plt.figure(1)
plt.title('Ag')
plt.ylabel(r'Total energy $[Ha]$')
plt.xlabel(r'Volume $[bohr]^3$')
plt.plot(xa,curve(xa),'-')
plt.plot(Volume,totalEnergy,'o')
plt.plot(min,curve(min),'o')
plt.annotate('optimal volume '+str(min[0]), xy=(min,curve(min)), xycoords='data' ,
              xytext=(min-10,curve(min)+0.002) ,  arrowprops=dict(arrowstyle="->"))
print 'optimal volume '+str(min[0])
 
plt.savefig('EV.png')
 
plt.show()

To execute this script, call:

python EV.py

Additional excercises

• Improve the accuracy of the fit by adding an additional data point.
• Calculate the lattice constant of the optimal volume.