High Pressure Thermoelasticity of Body-centered Cubic Tantalum
O. G\"{u}lseren$^{(1,2,3)}$ and R.E. Cohen$^{(3)}$
We have investigated the thermoelasticity of body-centered cubic (bcc) tantalum from first principles by using the linearized augmented plane wave (LAPW) and mixed--basis pseudopotential methods for pressures up to 400~GPa and temperatures up to 10000~K. Electronic excitation contributions to the free energy were included from the band structures, and phonon contributions were included using the particle-in-a-cell (PIC) model. The computed elastic constants agree well with available ultrasonic and diamond anvil cell data at low pressures, and shock data at high pressures. The shear modulus $c_{44}$ and the anisotropy change behavior with increasing pressure around 150~GPa because of an electronic topological transition. We find that the main contribution of temperature to the elastic constants is from the thermal expansivity. The PIC model in conjunction with fast self-consistent techniques is shown to be a tractable approach to studying thermoelasticity.
psfig.tex
% Psfig/TeX 
\def\PsfigVersion{1.9}
% dvips version
%
% All psfig/tex software, documentation, and related files
% in this distribution of psfig/tex are 
% Copyright 1987, 1988, 1991 Trevor J. Darrell
%
% Permission is granted for use and non-profit distribution of psfig/tex 
% providing that this notice is clearly maintained. The right to
% distribute any portion of psfig/tex for profit or as part of any commercial
% product is specifically reserved for the author(s) of that portion.
%
% *** Feel free to make local modifications of psfig as you wish,
% *** but DO NOT post any changed or modified versions of ``psfig''
% *** directly to the net. Send them to me and I'll try to incorporate
% *** them into future versions. If you want to take the psfig code 
% *** and make a new program (subject to the copyright above), distribute it, 
% *** (and maintain it) that's fine, just don't call it psfig.
%
% Bugs and improvements to trevor@media.mit.edu.
%
% Thanks to Greg Hager (GDH) and Ned Batchelder for their contributions
% to the original version of this project.
%
% Modified by J. Daniel Smith on 9 October 1990 to accept the
% %%BoundingBox: comment with or without a space after the colon.  Stole
% file reading code from Tom Rokicki's EPSF.TEX file (see below).
%
% More modifications by J. Daniel Smith on 29 March 1991 to allow the
% the included PostScript figure to be rotated.  The amount of
% rotation is specified by the "angle=" parameter of the \psfig command.
%
% Modified by Robert Russell on June 25, 1991 to allow users to specify
% .ps filenames which don't yet exist, provided they explicitly provide
% boundingbox information via the \psfig command. Note: This will only work
% if the "file=" parameter follows all four "bb???=" parameters in the
% command. This is due to the order in which psfig interprets these params.
%
%  3 Jul 1991	JDS	check if file already read in once
%  4 Sep 1991	JDS	fixed incorrect computation of rotated
%			bounding box
% 25 Sep 1991	GVR	expanded synopsis of \psfig
% 14 Oct 1991	JDS	\fbox code from LaTeX so \psdraft works with TeX
%			changed \typeout to \ps@typeout
% 17 Oct 1991	JDS	added \psscalefirst and \psrotatefirst
%

% From: gvr@cs.brown.edu (George V. Reilly)
%
% \psdraft	draws an outline box, but doesn't include the figure
%		in the DVI file.  Useful for previewing.
%
% \psfull	includes the figure in the DVI file (default).
%
% \psscalefirst width= or height= specifies the size of the figure
% 		before rotation.
% \psrotatefirst (default) width= or height= specifies the size of the
% 		 figure after rotation.  Asymetric figures will
% 		 appear to shrink.
%
% \psfigurepath#1	sets the path to search for the figure
%
% \psfig
% usage: \psfig{file=, figure=, height=, width=,
%			bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=, angle=, silent=}
%
%	"file" is the filename.  If no path name is specified and the
%		file is not found in the current directory,
%		it will be looked for in directory \psfigurepath.
%	"figure" is a synonym for "file".
%	By default, the width and height of the figure are taken from
%		the BoundingBox of the figure.
%	If "width" is specified, the figure is scaled so that it has
%		the specified width.  Its height changes proportionately.
%	If "height" is specified, the figure is scaled so that it has
%		the specified height.  Its width changes proportionately.
%	If both "width" and "height" are specified, the figure is scaled
%		anamorphically.
%	"bbllx", "bblly", "bburx", and "bbury" control the PostScript
%		BoundingBox.  If these four values are specified
%               *before* the "file" option, the PSFIG will not try to
%               open the PostScript file.
%	"rheight" and "rwidth" are the reserved height and width
%		of the figure, i.e., how big TeX actually thinks
%		the figure is.  They default to "width" and "height".
%	The "clip" option ensures that no portion of the figure will
%		appear outside its BoundingBox.  "clip=" is a switch and
%		takes no value, but the `=' must be present.
%	The "angle" option specifies the angle of rotation (degrees, ccw).
%	The "silent" option makes \psfig work silently.
%

% check to see if macros already loaded in (maybe some other file says
% "\input psfig") ...
\ifx\undefined\psfig\else\endinput\fi

%
% from a suggestion by eijkhout@csrd.uiuc.edu to allow
% loading as a style file. Changed to avoid problems
% with amstex per suggestion by jbence@math.ucla.edu

\let\LaTeXAtSign=\@
\let\@=\relax
\edef\psfigRestoreAt{\catcode`\@=\number\catcode`@\relax}
%\edef\psfigRestoreAt{\catcode`@=\number\catcode`@\relax}
\catcode`\@=11\relax
\newwrite\@unused
\def\ps@typeout#1{{\let\protect\string\immediate\write\@unused{#1}}}
\ps@typeout{psfig/tex \PsfigVersion}

%% Here's how you define your figure path.  Should be set up with null
%% default and a user useable definition.

\def\figurepath{./}
\def\psfigurepath#1{\edef\figurepath{#1}}

%
% @psdo control structure -- similar to Latex @for.
% I redefined these with different names so that psfig can
% be used with TeX as well as LaTeX, and so that it will not 
% be vunerable to future changes in LaTeX's internal
% control structure,
%
\def\@nnil{\@nil}
\def\@empty{}
\def\@psdonoop#1\@@#2#3{}
\def\@psdo#1:=#2\do#3{\edef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \expandafter\@psdoloop#2,\@nil,\@nil\@@#1{#3}\fi}
\def\@psdoloop#1,#2,#3\@@#4#5{\def#4{#1}\ifx #4\@nnil \else
       #5\def#4{#2}\ifx #4\@nnil \else#5\@ipsdoloop #3\@@#4{#5}\fi\fi}
\def\@ipsdoloop#1,#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@ipsdoloop\fi\@nextwhile#2\@@#3{#4}}
\def\@tpsdo#1:=#2\do#3{\xdef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \@tpsdoloop#2\@nil\@nil\@@#1{#3}\fi}
\def\@tpsdoloop#1#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@tpsdoloop\fi\@nextwhile#2\@@#3{#4}}
% 
% \fbox is defined in latex.tex; so if \fbox is undefined, assume that
% we are not in LaTeX.
% Perhaps this could be done better???
\ifx\undefined\fbox
% \fbox code from modified slightly from LaTeX
\newdimen\fboxrule
\newdimen\fboxsep
\newdimen\ps@tempdima
\newbox\ps@tempboxa
\fboxsep = 3pt
\fboxrule = .4pt
\long\def\fbox#1{\leavevmode\setbox\ps@tempboxa\hbox{#1}\ps@tempdima\fboxrule
    \advance\ps@tempdima \fboxsep \advance\ps@tempdima \dp\ps@tempboxa
   \hbox{\lower \ps@tempdima\hbox
  {\vbox{\hrule height \fboxrule
          \hbox{\vrule width \fboxrule \hskip\fboxsep
          \vbox{\vskip\fboxsep \box\ps@tempboxa\vskip\fboxsep}\hskip 
                 \fboxsep\vrule width \fboxrule}
                 \hrule height \fboxrule}}}}
\fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% file reading stuff from epsf.tex
%   EPSF.TEX macro file:
%   Written by Tomas Rokicki of Radical Eye Software, 29 Mar 1989.
%   Revised by Don Knuth, 3 Jan 1990.
%   Revised by Tomas Rokicki to accept bounding boxes with no
%      space after the colon, 18 Jul 1990.
%   Portions modified/removed for use in PSFIG package by
%      J. Daniel Smith, 9 October 1990.
%
\newread\ps@stream
\newif\ifnot@eof       % continue looking for the bounding box?
\newif\if@noisy        % report what you're making?
\newif\if@atend        % %%BoundingBox: has (at end) specification
\newif\if@psfile       % does this look like a PostScript file?
%
% PostScript files should start with `%!'
%
{\catcode`\%=12\global\gdef\epsf@start{%!}}
\def\epsf@PS{PS}
%
\def\epsf@getbb#1{%
%
%   The first thing we need to do is to open the
%   PostScript file, if possible.
%
\openin\ps@stream=#1
\ifeof\ps@stream\ps@typeout{Error, File #1 not found}\else
%
%   Okay, we got it. Now we'll scan lines until we find one that doesn't
%   start with %. We're looking for the bounding box comment.
%
   {\not@eoftrue \chardef\other=12
    \def\do##1{\catcode`##1=\other}\dospecials \catcode`\ =10
    \loop
       \if@psfile
	  \read\ps@stream to \epsf@fileline
       \else{
	  \obeyspaces
          \read\ps@stream to \epsf@tmp\global\let\epsf@fileline\epsf@tmp}
       \fi
       \ifeof\ps@stream\not@eoffalse\else
%
%   Check the first line for `%!'.  Issue a warning message if its not
%   there, since the file might not be a PostScript file.
%
       \if@psfile\else
       \expandafter\epsf@test\epsf@fileline:. \\%
       \fi
%
%   We check to see if the first character is a % sign;
%   if so, we look further and stop only if the line begins with
%   `%%BoundingBox:' and the `(atend)' specification was not found.
%   That is, the only way to stop is when the end of file is reached,
%   or a `%%BoundingBox: llx lly urx ury' line is found.
%
          \expandafter\epsf@aux\epsf@fileline:. \\%
       \fi
   \ifnot@eof\repeat
   }\closein\ps@stream\fi}%
%
% This tests if the file we are reading looks like a PostScript file.
%
\long\def\epsf@test#1#2#3:#4\\{\def\epsf@testit{#1#2}
			\ifx\epsf@testit\epsf@start\else
\ps@typeout{Warning! File does not start with `\epsf@start'.  It may not be a PostScript file.}
			\fi
			\@psfiletrue} % don't test after 1st line
%
%   We still need to define the tricky \epsf@aux macro. This requires
%   a couple of magic constants for comparison purposes.
%
{\catcode`\%=12\global\let\epsf@percent=%\global\def\epsf@bblit{%BoundingBox}}
%
%
%   So we're ready to check for `%BoundingBox:' and to grab the
%   values if they are found.  We continue searching if `(at end)'
%   was found after the `%BoundingBox:'.
%
\long\def\epsf@aux#1#2:#3\\{\ifx#1\epsf@percent
   \def\epsf@testit{#2}\ifx\epsf@testit\epsf@bblit
	\@atendfalse
        \epsf@atend #3 . \\%
	\if@atend	
	   \if@verbose{
		\ps@typeout{psfig: found `(atend)'; continuing search}
	   }\fi
        \else
        \epsf@grab #3 . . . \\%
        \not@eoffalse
        \global\no@bbfalse
        \fi
   \fi\fi}%
%
%   Here we grab the values and stuff them in the appropriate definitions.
%
\def\epsf@grab #1 #2 #3 #4 #5\\{%
   \global\def\epsf@llx{#1}\ifx\epsf@llx\empty
      \epsf@grab #2 #3 #4 #5 .\\\else
   \global\def\epsf@lly{#2}%
   \global\def\epsf@urx{#3}\global\def\epsf@ury{#4}\fi}%
%
% Determine if the stuff following the %%BoundingBox is `(atend)'
% J. Daniel Smith.  Copied from \epsf@grab above.
%
\def\epsf@atendlit{(atend)} 
\def\epsf@atend #1 #2 #3\\{%
   \def\epsf@tmp{#1}\ifx\epsf@tmp\empty
      \epsf@atend #2 #3 .\\\else
   \ifx\epsf@tmp\epsf@atendlit\@atendtrue\fi\fi}


% End of file reading stuff from epsf.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trigonometry stuff from "trig.tex"
\chardef\psletter = 11 % won't conflict with \begin{letter} now...
\chardef\other = 12

\newif \ifdebug %%% turn me on to see TeX hard at work ...
\newif\ifc@mpute %%% don't need to compute some values
\c@mputetrue % but assume that we do

\let\then = \relax
\def\r@dian{pt }
\let\r@dians = \r@dian
\let\dimensionless@nit = \r@dian
\let\dimensionless@nits = \dimensionless@nit
\def\internal@nit{sp }
\let\internal@nits = \internal@nit
\newif\ifstillc@nverging
\def \Mess@ge #1{\ifdebug \then \message {#1} \fi}

{ %%% Things that need abnormal catcodes %%%
	\catcode `\@ = \psletter
	\gdef \nodimen {\expandafter \n@dimen \the \dimen}
	\gdef \term #1 #2 #3%
	       {\edef \t@ {\the #1}%%% freeze parameter 1 (count, by value)
		\edef \t@@ {\expandafter \n@dimen \the #2\r@dian}%
				   %%% freeze parameter 2 (dimen, by value)
		\t@rm {\t@} {\t@@} {#3}%
	       }
	\gdef \t@rm #1 #2 #3%
	       {{%
		\count 0 = 0
		\dimen 0 = 1 \dimensionless@nit
		\dimen 2 = #2\relax
		\Mess@ge {Calculating term #1 of \nodimen 2}%
		\loop
		\ifnum	\count 0 < #1
		\then	\advance \count 0 by 1
			\Mess@ge {Iteration \the \count 0 \space}%
			\Multiply \dimen 0 by {\dimen 2}%
			\Mess@ge {After multiplication, term = \nodimen 0}%
			\Divide \dimen 0 by {\count 0}%
			\Mess@ge {After division, term = \nodimen 0}%
		\repeat
		\Mess@ge {Final value for term #1 of 
				\nodimen 2 \space is \nodimen 0}%
		\xdef \Term {#3 = \nodimen 0 \r@dians}%
		\aftergroup \Term
	       }}
	\catcode `\p = \other
	\catcode `\t = \other
	\gdef \n@dimen #1pt{#1} %%% throw away the ``pt''
}

\def \Divide #1by #2{\divide #1 by #2} %%% just a synonym

\def \Multiply #1by #2%%% allows division of a dimen by a dimen
       {{%%% should really freeze parameter 2 (dimen, passed by value)
	\count 0 = #1\relax
	\count 2 = #2\relax
	\count 4 = 65536
	\Mess@ge {Before scaling, count 0 = \the \count 0 \space and
			count 2 = \the \count 2}%
	\ifnum	\count 0 > 32767 %%% do our best to avoid overflow
	\then	\divide \count 0 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 0 < -32767
		\then	\divide \count 0 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\ifnum	\count 2 > 32767 %%% while retaining reasonable accuracy
	\then	\divide \count 2 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 2 < -32767
		\then	\divide \count 2 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\multiply \count 0 by \count 2
	\divide \count 0 by \count 4
	\xdef \product {#1 = \the \count 0 \internal@nits}%
	\aftergroup \product
       }}

\def\r@duce{\ifdim\dimen0 > 90\r@dian \then   % sin(x+90) = sin(180-x)
		\multiply\dimen0 by -1
		\advance\dimen0 by 180\r@dian
		\r@duce
	    \else \ifdim\dimen0 < -90\r@dian \then  % sin(-x) = sin(360+x)
		\advance\dimen0 by 360\r@dian
		\r@duce
		\fi
	    \fi}

\def\Sine#1%
       {{%
	\dimen 0 = #1 \r@dian
	\r@duce
	\ifdim\dimen0 = -90\r@dian \then
	   \dimen4 = -1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 90\r@dian \then
	   \dimen4 = 1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 0\r@dian \then
	   \dimen4 = 0\r@dian
	   \c@mputefalse
	\fi
%
	\ifc@mpute \then
        	% convert degrees to radians
		\divide\dimen0 by 180
		\dimen0=3.141592654\dimen0
%
		\dimen 2 = 3.1415926535897963\r@dian %%% a well-known constant
		\divide\dimen 2 by 2 %%% we only deal with -pi/2 : pi/2
		\Mess@ge {Sin: calculating Sin of \nodimen 0}%
		\count 0 = 1 %%% see power-series expansion for sine
		\dimen 2 = 1 \r@dian %%% ditto
		\dimen 4 = 0 \r@dian %%% ditto
		\loop
			\ifnum	\dimen 2 = 0 %%% then we've done
			\then	\stillc@nvergingfalse 
			\else	\stillc@nvergingtrue
			\fi
			\ifstillc@nverging %%% then calculate next term
			\then	\term {\count 0} {\dimen 0} {\dimen 2}%
				\advance \count 0 by 2
				\count 2 = \count 0
				\divide \count 2 by 2
				\ifodd	\count 2 %%% signs alternate
				\then	\advance \dimen 4 by \dimen 2
				\else	\advance \dimen 4 by -\dimen 2
				\fi
		\repeat
	\fi		
			\xdef \sine {\nodimen 4}%
       }}

% Now the Cosine can be calculated easily by calling \Sine
\def\Cosine#1{\ifx\sine\UnDefined\edef\Savesine{\relax}\else
		             \edef\Savesine{\sine}\fi
	{\dimen0=#1\r@dian\advance\dimen0 by 90\r@dian
	 \Sine{\nodimen 0}
	 \xdef\cosine{\sine}
	 \xdef\sine{\Savesine}}}	      
% end of trig stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\psdraft{
	\def\@psdraft{0}
	%\ps@typeout{draft level now is \@psdraft \space . }
}
\def\psfull{
	\def\@psdraft{100}
	%\ps@typeout{draft level now is \@psdraft \space . }
}

\psfull

\newif\if@scalefirst
\def\psscalefirst{\@scalefirsttrue}
\def\psrotatefirst{\@scalefirstfalse}
\psrotatefirst

\newif\if@draftbox
\def\psnodraftbox{
	\@draftboxfalse
}
\def\psdraftbox{
	\@draftboxtrue
}
\@draftboxtrue

\newif\if@prologfile
\newif\if@postlogfile
\def\pssilent{
	\@noisyfalse
}
\def\psnoisy{
	\@noisytrue
}
\psnoisy
%%% These are for the option list.
%%% A specification of the form a = b maps to calling \@p@@sa{b}
\newif\if@bbllx
\newif\if@bblly
\newif\if@bburx
\newif\if@bbury
\newif\if@height
\newif\if@width
\newif\if@rheight
\newif\if@rwidth
\newif\if@angle
\newif\if@clip
\newif\if@verbose
\def\@p@@sclip#1{\@cliptrue}


\newif\if@decmpr

%%% GDH 7/26/87 -- changed so that it first looks in the local directory,
%%% then in a specified global directory for the ps file.
%%% RPR 6/25/91 -- changed so that it defaults to user-supplied name if
%%% boundingbox info is specified, assuming graphic will be created by
%%% print time.
%%% TJD 10/19/91 -- added bbfile vs. file distinction, and @decmpr flag

\def\@p@@sfigure#1{\def\@p@sfile{null}\def\@p@sbbfile{null}
	        \openin1=#1.bb
		\ifeof1\closein1
	        	\openin1=\figurepath#1.bb
			\ifeof1\closein1
			        \openin1=#1
				\ifeof1\closein1%
				       \openin1=\figurepath#1
					\ifeof1
					   \ps@typeout{Error, File #1 not found}
						\if@bbllx\if@bblly
				   		\if@bburx\if@bbury
			      				\def\@p@sfile{#1}%
			      				\def\@p@sbbfile{#1}%
							\@decmprfalse
				  	   	\fi\fi\fi\fi
					\else\closein1
				    		\def\@p@sfile{\figurepath#1}%
				    		\def\@p@sbbfile{\figurepath#1}%
						\@decmprfalse
	                       		\fi%
			 	\else\closein1%
					\def\@p@sfile{#1}
					\def\@p@sbbfile{#1}
					\@decmprfalse
			 	\fi
			\else
				\def\@p@sfile{\figurepath#1}
				\def\@p@sbbfile{\figurepath#1.bb}
				\@decmprtrue
			\fi
		\else
			\def\@p@sfile{#1}
			\def\@p@sbbfile{#1.bb}
			\@decmprtrue
		\fi}

\def\@p@@sfile#1{\@p@@sfigure{#1}}

\def\@p@@sbbllx#1{
		%\ps@typeout{bbllx is #1}
		\@bbllxtrue
		\dimen100=#1
		\edef\@p@sbbllx{\number\dimen100}
}
\def\@p@@sbblly#1{
		%\ps@typeout{bblly is #1}
		\@bbllytrue
		\dimen100=#1
		\edef\@p@sbblly{\number\dimen100}
}
\def\@p@@sbburx#1{
		%\ps@typeout{bburx is #1}
		\@bburxtrue
		\dimen100=#1
		\edef\@p@sbburx{\number\dimen100}
}
\def\@p@@sbbury#1{
		%\ps@typeout{bbury is #1}
		\@bburytrue
		\dimen100=#1
		\edef\@p@sbbury{\number\dimen100}
}
\def\@p@@sheight#1{
		\@heighttrue
		\dimen100=#1
   		\edef\@p@sheight{\number\dimen100}
		%\ps@typeout{Height is \@p@sheight}
}
\def\@p@@swidth#1{
		%\ps@typeout{Width is #1}
		\@widthtrue
		\dimen100=#1
		\edef\@p@swidth{\number\dimen100}
}
\def\@p@@srheight#1{
		%\ps@typeout{Reserved height is #1}
		\@rheighttrue
		\dimen100=#1
		\edef\@p@srheight{\number\dimen100}
}
\def\@p@@srwidth#1{
		%\ps@typeout{Reserved width is #1}
		\@rwidthtrue
		\dimen100=#1
		\edef\@p@srwidth{\number\dimen100}
}
\def\@p@@sangle#1{
		%\ps@typeout{Rotation is #1}
		\@angletrue
%		\dimen100=#1
		\edef\@p@sangle{#1} %\number\dimen100}
}
\def\@p@@ssilent#1{ 
		\@verbosefalse
}
\def\@p@@sprolog#1{\@prologfiletrue\def\@prologfileval{#1}}
\def\@p@@spostlog#1{\@postlogfiletrue\def\@postlogfileval{#1}}
\def\@cs@name#1{\csname #1\endcsname}
\def\@setparms#1=#2,{\@cs@name{@p@@s#1}{#2}}
%
% initialize the defaults (size the size of the figure)
%
\def\ps@init@parms{
		\@bbllxfalse \@bbllyfalse
		\@bburxfalse \@bburyfalse
		\@heightfalse \@widthfalse
		\@rheightfalse \@rwidthfalse
		\def\@p@sbbllx{}\def\@p@sbblly{}
		\def\@p@sbburx{}\def\@p@sbbury{}
		\def\@p@sheight{}\def\@p@swidth{}
		\def\@p@srheight{}\def\@p@srwidth{}
		\def\@p@sangle{0}
		\def\@p@sfile{} \def\@p@sbbfile{}
		\def\@p@scost{10}
		\def\@sc{}
		\@prologfilefalse
		\@postlogfilefalse
		\@clipfalse
		\if@noisy
			\@verbosetrue
		\else
			\@verbosefalse
		\fi
}
%
% Go through the options setting things up.
%
\def\parse@ps@parms#1{
	 	\@psdo\@psfiga:=#1\do
		   {\expandafter\@setparms\@psfiga,}}
%
% Compute bb height and width
%
\newif\ifno@bb
\def\bb@missing{
	\if@verbose{
		\ps@typeout{psfig: searching \@p@sbbfile \space  for bounding box}
	}\fi
	\no@bbtrue
	\epsf@getbb{\@p@sbbfile}
        \ifno@bb \else \bb@cull\epsf@llx\epsf@lly\epsf@urx\epsf@ury\fi
}	
\def\bb@cull#1#2#3#4{
	\dimen100=#1 bp\edef\@p@sbbllx{\number\dimen100}
	\dimen100=#2 bp\edef\@p@sbblly{\number\dimen100}
	\dimen100=#3 bp\edef\@p@sbburx{\number\dimen100}
	\dimen100=#4 bp\edef\@p@sbbury{\number\dimen100}
	\no@bbfalse
}
% rotate point (#1,#2) about (0,0).
% The sine and cosine of the angle are already stored in \sine and
% \cosine.  The result is placed in (\p@intvaluex, \p@intvaluey).
\newdimen\p@intvaluex
\newdimen\p@intvaluey
\def\rotate@#1#2{{\dimen0=#1 sp\dimen1=#2 sp
%            	calculate x' = x \cos\theta - y \sin\theta
		  \global\p@intvaluex=\cosine\dimen0
		  \dimen3=\sine\dimen1
		  \global\advance\p@intvaluex by -\dimen3
% 		calculate y' = x \sin\theta + y \cos\theta
		  \global\p@intvaluey=\sine\dimen0
		  \dimen3=\cosine\dimen1
		  \global\advance\p@intvaluey by \dimen3
		  }}
\def\compute@bb{
		\no@bbfalse
		\if@bbllx \else \no@bbtrue \fi
		\if@bblly \else \no@bbtrue \fi
		\if@bburx \else \no@bbtrue \fi
		\if@bbury \else \no@bbtrue \fi
		\ifno@bb \bb@missing \fi
		\ifno@bb \ps@typeout{FATAL ERROR: no bb supplied or found}
			\no-bb-error
		\fi
		%
%\ps@typeout{BB: \@p@sbbllx, \@p@sbblly, \@p@sbburx, \@p@sbbury} 
%
% store height/width of original (unrotated) bounding box
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\ps@bbw{\number\count203}
		\edef\ps@bbh{\number\count204}
		%\ps@typeout{ psbbh = \ps@bbh, psbbw = \ps@bbw }
		\if@angle 
			\Sine{\@p@sangle}\Cosine{\@p@sangle}
	        	{\dimen100=\maxdimen\xdef\r@p@sbbllx{\number\dimen100}
					    \xdef\r@p@sbblly{\number\dimen100}
			                    \xdef\r@p@sbburx{-\number\dimen100}
					    \xdef\r@p@sbbury{-\number\dimen100}}
%
% Need to rotate all four points and take the X-Y extremes of the new
% points as the new bounding box.
                        \def\minmaxtest{
			   \ifnum\number\p@intvaluex<\r@p@sbbllx
			      \xdef\r@p@sbbllx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluex>\r@p@sbburx
			      \xdef\r@p@sbburx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluey<\r@p@sbblly
			      \xdef\r@p@sbblly{\number\p@intvaluey}\fi
			   \ifnum\number\p@intvaluey>\r@p@sbbury
			      \xdef\r@p@sbbury{\number\p@intvaluey}\fi
			   }
%			lower left
			\rotate@{\@p@sbbllx}{\@p@sbblly}
			\minmaxtest
%			upper left
			\rotate@{\@p@sbbllx}{\@p@sbbury}
			\minmaxtest
%			lower right
			\rotate@{\@p@sbburx}{\@p@sbblly}
			\minmaxtest
%			upper right
			\rotate@{\@p@sbburx}{\@p@sbbury}
			\minmaxtest
			\edef\@p@sbbllx{\r@p@sbbllx}\edef\@p@sbblly{\r@p@sbblly}
			\edef\@p@sbburx{\r@p@sbburx}\edef\@p@sbbury{\r@p@sbbury}
%\ps@typeout{rotated BB: \r@p@sbbllx, \r@p@sbblly, \r@p@sbburx, \r@p@sbbury}
		\fi
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\@bbw{\number\count203}
		\edef\@bbh{\number\count204}
		%\ps@typeout{ bbh = \@bbh, bbw = \@bbw }
}
%
% \in@hundreds performs #1 * (#2 / #3) correct to the hundreds,
%	then leaves the result in @result
%
\def\in@hundreds#1#2#3{\count240=#2 \count241=#3
		     \count100=\count240	% 100 is first digit #2/#3
		     \divide\count100 by \count241
		     \count101=\count100
		     \multiply\count101 by \count241
		     \advance\count240 by -\count101
		     \multiply\count240 by 10
		     \count101=\count240	%101 is second digit of #2/#3
		     \divide\count101 by \count241
		     \count102=\count101
		     \multiply\count102 by \count241
		     \advance\count240 by -\count102
		     \multiply\count240 by 10
		     \count102=\count240	% 102 is the third digit
		     \divide\count102 by \count241
		     \count200=#1\count205=0
		     \count201=\count200
			\multiply\count201 by \count100
		 	\advance\count205 by \count201
		     \count201=\count200
			\divide\count201 by 10
			\multiply\count201 by \count101
			\advance\count205 by \count201
			%
		     \count201=\count200
			\divide\count201 by 100
			\multiply\count201 by \count102
			\advance\count205 by \count201
			%
		     \edef\@result{\number\count205}
}
\def\compute@wfromh{
		% computing : width = height * (bbw / bbh)
		\in@hundreds{\@p@sheight}{\@bbw}{\@bbh}
		%\ps@typeout{ \@p@sheight * \@bbw / \@bbh, = \@result }
		\edef\@p@swidth{\@result}
		%\ps@typeout{w from h: width is \@p@swidth}
}
\def\compute@hfromw{
		% computing : height = width * (bbh / bbw)
	        \in@hundreds{\@p@swidth}{\@bbh}{\@bbw}
		%\ps@typeout{ \@p@swidth * \@bbh / \@bbw = \@result }
		\edef\@p@sheight{\@result}
		%\ps@typeout{h from w : height is \@p@sheight}
}
\def\compute@handw{
		\if@height 
			\if@width
			\else
				\compute@wfromh
			\fi
		\else 
			\if@width
				\compute@hfromw
			\else
				\edef\@p@sheight{\@bbh}
				\edef\@p@swidth{\@bbw}
			\fi
		\fi
}
\def\compute@resv{
		\if@rheight \else \edef\@p@srheight{\@p@sheight} \fi
		\if@rwidth \else \edef\@p@srwidth{\@p@swidth} \fi
		%\ps@typeout{rheight = \@p@srheight, rwidth = \@p@srwidth}
}
%		
% Compute any missing values
\def\compute@sizes{
	\compute@bb
	\if@scalefirst\if@angle
% at this point the bounding box has been adjsuted correctly for
% rotation.  PSFIG does all of its scaling using \@bbh and \@bbw.  If
% a width= or height= was specified along with \psscalefirst, then the
% width=/height= value needs to be adjusted to match the new (rotated)
% bounding box size (specifed in \@bbw and \@bbh).
%    \ps@bbw       width=
%    -------  =  ---------- 
%    \@bbw       new width=
% so `new width=' = (width= * \@bbw) / \ps@bbw; where \ps@bbw is the
% width of the original (unrotated) bounding box.
	\if@width
	   \in@hundreds{\@p@swidth}{\@bbw}{\ps@bbw}
	   \edef\@p@swidth{\@result}
	\fi
	\if@height
	   \in@hundreds{\@p@sheight}{\@bbh}{\ps@bbh}
	   \edef\@p@sheight{\@result}
	\fi
	\fi\fi
	\compute@handw
	\compute@resv}

%
% \psfig
% usage : \psfig{file=, height=, width=, bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=}
%
% "clip=" is a switch and takes no value, but the `=' must be present.
\def\psfig#1{\vbox {
	% do a zero width hard space so that a single
	% \psfig in a centering enviornment will behave nicely
	%{\setbox0=\hbox{\ }\ \hskip-\wd0}
	%
	\ps@init@parms
	\parse@ps@parms{#1}
	\compute@sizes
	%
	\ifnum\@p@scost<\@psdraft{
		%
		\special{ps::[begin] 	\@p@swidth \space \@p@sheight \space
				\@p@sbbllx \space \@p@sbblly \space
				\@p@sbburx \space \@p@sbbury \space
				startTexFig \space }
		\if@angle
			\special {ps:: \@p@sangle \space rotate \space} 
		\fi
		\if@clip{
			\if@verbose{
				\ps@typeout{(clip)}
			}\fi
			\special{ps:: doclip \space }
		}\fi
		\if@prologfile
		    \special{ps: plotfile \@prologfileval \space } \fi
		\if@decmpr{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile.Z \space }
			}\fi
			\special{ps: plotfile "`zcat \@p@sfile.Z" \space }
		}\else{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile \space }
			}\fi
			\special{ps: plotfile \@p@sfile \space }
		}\fi
		\if@postlogfile
		    \special{ps: plotfile \@postlogfileval \space } \fi
		\special{ps::[end] endTexFig \space }
		% Create the vbox to reserve the space for the figure.
		\vbox to \@p@srheight sp{
		% 1/92 TJD Changed from "true sp" to "sp" for magnification.
			\hbox to \@p@srwidth sp{
				\hss
			}
		\vss
		}
	}\else{
		% draft figure, just reserve the space and print the
		% path name.
		\if@draftbox{		
			% Verbose draft: print file name in box
			\hbox{\frame{\vbox to \@p@srheight sp{
			\vss
			\hbox to \@p@srwidth sp{ \hss \@p@sfile \hss }
			\vss
			}}}
		}\else{
			% Non-verbose draft
			\vbox to \@p@srheight sp{
			\vss
			\hbox to \@p@srwidth sp{\hss}
			\vss
			}
		}\fi	



	}\fi
}}
\psfigRestoreAt
\let\@=\LaTeXAtSign

taecons2col.tex
\documentstyle[twocolumn,prb,aps]{revtex}
%\documentstyle[preprint,prb,aps]{revtex}

\begin{document}
\draft
\input{psfig}
\preprint{ }
% FOR TWO COLUMN ACTIVATE THE LINE BELOW 
\twocolumn[\hsize\textwidth\columnwidth\hsize\csname @twocolumnfalse\endcsname

\title{High Pressure Thermoelasticity of Body-centered Cubic Tantalum}

\author{O. G\"{u}lseren$^{(1,2,3)}$ and R.E. Cohen$^{(3)}$}
\address{$^{(1)}$ NIST Center for Neutron Research,
National Institute of Standards and Technology, Gaithersburg, MD 20899}
\address{$^{(2)}$ Department of Materials Science and Engineering,
University of Pennsylvania, Philadelphia, PA 19104}
\address{$^{(3)}$
Geophysical Laboratory and Center for High Pressure Research,
Carnegie Institution of Washington,
5251 Broad Branch Road, NW, Washington, DC 20015.}

\date{\today}
\maketitle

\begin{abstract}
We have investigated the thermoelasticity of body-centered cubic (bcc)
tantalum from first principles by using the linearized augmented plane
wave (LAPW) and mixed--basis pseudopotential methods for pressures up
to 400~GPa and temperatures up to 10000~K.
Electronic excitation contributions to the free energy were included
from the band structures, and phonon contributions were included using
the particle-in-a-cell (PIC) model. The computed elastic constants agree
well with available ultrasonic and diamond anvil cell data at low
pressures, and shock data at high pressures. The shear modulus $c_{44}$
and the anisotropy change behavior with increasing pressure around
150~GPa because of an electronic topological transition. We find that
the main contribution of temperature to the elastic constants is from
the thermal expansivity. The PIC model in conjunction with fast
self-consistent techniques is shown to be a tractable approach to
studying thermoelasticity.
\end{abstract}

%\vspace{2cm}
\pacs{PACS numbers: 62.20.Dc, 62.20.-x, 71.20.Be, 46.25.Hf, 46.25.-y}

%62.20.Dc Elasticity, elastic constants
%62.20.-x Mechanical properties of solids 
%71.20.Be Transition metals and alloys
%46.25.Hf Thermoelasticity and electromagnetic elasticity
%46.25.-y Static elasticity
%
%05.70.Ce Thermodynamic functions and equations of state
%64.30.+t Equations of state of specific substances
%65.50.+m Thermodynamic properties and entropy
%65.70.+y Thermal expansion and density changes; thermomechanical effects

]

Single crystal elastic constants of solids at high pressures and
temperatures are essential in order to predict and understand material
response, strength, mechanical stability, and phase transitions.
We have studied the high pressure and temperature elastic constants of 
body-centered cubic (bcc) tantalum, a group V transition metal, from
{\sf first principles}. Because of its high structural mechanical, thermal
and chemical stability, Ta is a useful high pressure standard.\cite{cynn}
Ta has a very high melting temperature, 3269~K at ambient pressure.
Bcc Ta is stable to 174~GPa, according to diamond-anvil-cell
experiments.\cite{cynn} Shock compression experiments~\cite{hug1} show no
transition other than melting (at around 300~GPa). Its stability makes Ta
an ideal material for understanding the generic behavior of transition
metals under compression, without the complication of phase transitions.
Recently, its static properties were studied by full-potential LMTO
calculations~\cite{soderlind} and the thermal equation of state was
reported.\cite{tacohen}

The three elastic constants, $c_{11}$, $c_{12}$ and $c_{44}$, completely
describe the elastic behaviour of a cubic crystal. A more convenient set
for computations are $c_{44}$ and two linear combinations, $K$ and $c_s$.
The bulk modulus, 
\begin{equation}
K=(c_{11}+2c_{12})/3,
\end{equation}
is the resistance to deformation by an uniform hydrostatic pressure;
the shear constant,
\begin{equation}
c_s=(c_{11}-c_{12})/2,
\end{equation}
is the resistance to shear deformation across the (110) plane in the
$[1\bar{1}0]$ direction, and $c_{44}$ is the resistance to shear
deformation across the (100) plane in the $[010]$ direction.
The bulk modulus $K$ was determined from the equation of
state,\cite{tacohen} using the Vinet equation.\cite{vinet,ron1}
We obtained the shear moduli by straining the bcc lattice at fixed volumes
using volume conserving tetragonal and orthorhombic strains for $c_s$ and
$c_{44}$ respectively, and computing the free energy as a function of
strain. $c_s$ was obtained by applying the following isochoric strain
\begin{equation}
\epsilon = \left( \begin{array}{ccc}
\delta & 0 & 0 \\
0 & \delta & 0 \\
0 & 0 & (1+\delta)^{-2}-1
\end{array}
\right) ,
\end{equation}
where $\delta$ is the magnitude of the strain.
Then the strain energy is
\begin{equation}
F(\delta) = F(0) + 6 c_s V \delta^2 + O(\delta^3),
\end{equation}
where $F(0)$ is the free energy of the unstrained system and $V$ is its
volume. Similarly, $c_{44}$ was calculated from the following strain
\begin{equation}
\epsilon = \left( \begin{array}{ccc}
0 & \delta & 0 \\
\delta & 0 & 0 \\
0 & 0 & \delta^2/(1-\delta^2)
\end{array}
\right)
\end{equation}
with the corresponding strain energy
\begin{equation}
F(\delta) = F(0) + 2 c_{44} V \delta^2 + O(\delta^4).
\end{equation}
The quadratic coefficients of strain energy gives the elastic
constants. First order terms due to the initial stress
(hydrostatic pressure)~\cite{barron} were eliminated by applying isochoric
strains. Then, the elastic constants $c_{11}$ and $c_{12}$ were obtained from
$c_s$ and $K$.

We assume that the Helmholtz free energy of the system can be separated
as~\cite{fecons,wass,iron}:
\begin{equation}
F(V,T) = E_{0}(V) + F_{el}(V,T) + F_{vib}(V,T)
\end{equation}
where $E_{0}(V)$ is the static zero temperature energy,
$F_{el}(V,T)$ is the electronic contribution, and
$F_{vib}(V,T)$ is the vibrational contribution to the free energy.
Our computational procedure is based on density functional theory (DFT)
generalized to finite temperatures by the Mermin theorem.\cite{temdft}
The charge density is temperature dependent through occupation numbers
according to the Fermi-Dirac distribution, giving the electronic
entropy from
\begin{equation}
S_{el}=\sum{f_i ln f_i + (1-f_i) ln (1-f_i)}
\end{equation}
where $f_i=f_i(E-E_F,T))$ is the Fermi occupation at $T$ for each
state $i$. The variations of $f_i$ with temperature were included from
the self-consistent band structures calculated at an electronic
temperature of 2000~K varying according to the Fermi-Dirac distribution.

The electronic excitations, both the static energy and the electronic
contribution to the free energy, were computed by using the full potential
linearized augmented plane wave (LAPW) method.\cite{weikrakauer,singh}
The $5p$,$4f$,$5d$ and $6s$ states were treated as band states, and the
deeper states were treated as soft core electrons. The generalized gradient
approximation (GGA)\cite{gga-pbe} was used for the exchange-correlation
potential. The convergence of strain energies with respect to the Brillouin
zone integration was carefully checked by repeating the calculations for
16x16x16 and 24x24x24 meshes at V=16.82~\AA$^3${} and we found at most 2~GPa
(3 \%) difference both for $c_{s}$ and $c_{44}$. Hence, we used 16x16x16
special k-points meshes~\cite{monpack} in the full Brillouin zone giving
344 and 612 k-points within IBZ of tetragonal and orthorhombic lattice
respectively. The convergence parameter $R K_{max}$ was 9 giving about
1800 planewaves and 200 basis functions per atom at zero pressure.

The vibrational free-energy was obtained within the particle-in-a-cell
model (PIC)~\cite{cellmodel} by using an accurate pseudopotential
mixed-basis total energy method~\cite{mix-oguz} which is computationally
more efficient than the LAPW calculations. In PIC, an atom is displaced in
its Wigner-Seitz cell in the potential field of all the other atoms fixed
at their equilibrium positions, i.e. the ideal, static lattice except for
the wanderer atom. The partition function, and hence the free energy is
calculated from this potential energy surface via an integral over the
position of a single atom inside the Wigner-Seitz cell. The PIC model is
essentially an anharmonic Einstein model, and the 3N dimensional partition
function is reduced to a simple 3D integral.\cite{tacohen,wass}
The advantage of the cell model over lattice dynamics based on the
quasiharmonic approximation is that anharmonic contributions from the
potential-energy of the system have been included exactly without a
perturbation expansion. On the other hand, since we used the classical
partition function, and the interatomic correlations between the motions of
different atoms is ignored, it is only valid at high temperatures above the
Debye temperature (245 K in Ta.\cite{grigoriev}) Since the vacancy formation
energy is very high in Ta,\cite{sonali} spontaneous formation of defects is
only important after the melting temperature.

\begin{figure}
\centerline{\psfig{figure=econst0k.ps,angle=0,width=82mm}}
%\vspace*{1mm}
\caption{\small Static elastic constants of Ta as a function of pressure.
Open squares are ultrasonic experimental data of Katahara 
{\it et al}\protect\cite{katahara1,katahara2}, and the dotted lines
show the initial slopes. Open symbols are SAX data of Cynn and
Yoo\protect\cite{cynnecons}.
The anisotropy ratio is shown in the inset.}
\label{fig:econst0K}
\end{figure}

For the PIC computations, a supercell with 54 atoms was used. The
pseudopotential mixed-basis calculations were carried out on this 54 atoms
supercell using LDA~\cite{ldazun} for exchange-correlations effects
and 2x2x2 k-point mesh resulting 4 special $\vec{k}$ points for BZ
integrations. A semi-relativistic, nonlocal and norm-conserving
Troullier-Martins~\cite{tromartins} pseudopotential (with associated
pseudo-atomic orbitals) with non-linear core corrections~\cite{nonlinearcore}
was used for the Ta atoms as described in detail in our previous study
of thermal equation of state of Ta.\cite{tacohen} After checking the
energy convergence, 550~eV and 60~eV are used as planewave energy cutoffs
for the expansion of the pseudo-atomic orbitals as well as FFT grid and
low-energy plane waves for additional degrees of freedom in basis set
respectively. The canonical partition function was computed from the
potential energy surface as a function of displacements of wanderer atom
along special symmetry directions.\cite{wass,spdir} We used 2 and 4 special
directions for tetragonal and orthorhombic distortions, respectively,
which integrates exactly up to $l=6$ lattice harmonics.\cite{spdir} The
potential energy was calculated at 4-6 different displacements along each
of these special directions, and was fit to an even polynomials up to order
six. Details of the all computational parameters were described
previously.\cite{tacohen}

The static elastic constants as functions of pressure are presented in
Fig.~\ref{fig:econst0K} and Table~\ref{table:econs}. The zero pressure
values and initial slopes are
in good agreement with the ultrasonic experimental data of Katahara
{\it et al.}\cite{katahara1,katahara2} Similarly, comparison with recent SAX
(stress/angle-resolved x-ray diffraction) experimental data~\cite{cynnecons}
up to 105~GPa shows good agreement for $c_{11}$ and $c_{44}$. Likewise,
$c_{12}$ agrees well at low pressures, but deviates with increasing pressure.
This may be due to the assumed isostress condition for experimental
data analysis for all pressures, or due to the large uncertainty on measured
deviatoric stress at high pressures. Note that, the initial slope of
ultrasonic data agrees very well with our calculated $c_{12}$.
The anisotropy ratio, $A=c_{44}/c_{s}$ (inset Fig.~\ref{fig:econst0K})
first decreases from 1.57 to 0.9 with increasing pressure, and then its
slope reverses and it increases with increasing pressure. This is due to
the changes in $c_{44}$. The change in behavior of $c_{44}$ around 150 GPa is
due to the electronic transition evident in the equation of
state.\cite{tacohen} This indicates that elastic constants can be much more
sensitive to changes in the Fermi surface than the equation of state, where
the electronic transition was not apparent without examining small residuals
in the equation of state fit.

The elastic constants of Ta as functions of temperature at various
pressures are presented in Fig.~\ref{fig:econsther} and
Table~\ref{table:econsther}.
In order to compare with experimental data, the computed isothermal elastic
constants ($c_{ij}^{T}$) are converted to adiabatic constants ($c_{ij}^{S}$)
according to~\cite{davies}
\begin{equation}
c_{ij}^{S} =  c_{ij}^{T} + \frac{T}{\rho C_V} \lambda_i \lambda_j
\label{eq:isoadb}
\end{equation}
where $\lambda_i = \sum_k\alpha_k c_{ik}^{T}$, $\alpha_k$ is the linear
thermal expansion tensor, $C_V$ is the specific heat and $\rho$ is the
density. For cubic crystals, eq.~\ref{eq:isoadb} simplifies to
\begin{eqnarray}
c_{11}^{S} & = & c_{11}^{T} + \Delta \\
c_{12}^{S} & = & c_{12}^{T} + \Delta
\end{eqnarray}
where
\begin{equation}
\Delta = T (\alpha K_T)^2/(\rho C_V) = \rho C_V T \gamma^2 = 
T \alpha K_T \gamma
\end{equation}

\begin{table}[h]
\begin{tabular}{r r r r r r r}
%\hline
V (\AA$^3$) & Pressure & K & $c_{44}$ & $c_s$ & $c_{11}$ & $c_{12}$ \\
\hline
18.39 &  -0.76 &  187.89 &  66.30 &  42.12 &  244.05 &  159.82 \\
17.66 &   7.56 &  220.35 &  75.13 &  51.41 &  288.90 &  186.08 \\
16.82 &  19.35 &  263.82 &  82.73 &  62.89 &  347.68 &  221.90 \\
15.22 &  50.88 &  371.07 & 101.15 &  91.30 &  492.79 &  310.20 \\
13.01 & 125.70 &  598.02 & 129.18 & 143.78 &  789.73 &  502.17 \\
12.43 & 160.63 &  696.44 & 156.64 &        &         &         \\
11.67 & 202.01 &  808.93 & 194.45 &        &         &         \\
11.03 & 250.90 &  937.44 & 253.95 & 241.48 & 1259.41 &  776.46 \\
 9.83 & 376.54 & 1251.81 & 417.35 & 335.65 & 1699.34 & 1028.04 \\
 9.26 & 456.48 & 1443.14 & 538.35 & 391.79 & 1965.52 & 1181.94 \\
%\hline
\end{tabular}
\vspace*{0.1cm}
\caption{The static elastic constants for bcc Tantalum. All elastic
constants as well as pressure units are GPa.}
\label{table:econs}
\end{table}
\begin{figure}
\centerline{\psfig{figure=econstadb.ps,angle=0,width=80mm}}
%\vspace*{3mm}
\caption{\small The elastic constants of bcc Ta as a function of
temperature at different pressure from 0~GPa (lowest curve)
to 400~GPa (uppermost curve) with 50~GPa interval.
a) shear modulus $c_{s}$, b) shear modulus
$c_{44}$ and c) adiabatic bulk modulus $K_S$.
Dotted lines are the experimental data from Walker
et al.~\protect\cite{walker}.}
\label{fig:econsther}
\end{figure}
with $\alpha$ is the thermal expansion coefficient, $\gamma$ is the
Gr\"{u}neisen parameter, and $K_T$ is the isothermal bulk modulus. The
thermodynamic parameters were computed self-consistently from the thermal
equation of state.\cite{tacohen} The correction is zero for $c_{44}$ and
$c_{s}$. $\Delta$ increases with temperature but decreases with pressure;
at 3000~K it decreases from 5~\% to 1~\% for pressures 50~GPa to 400~GPa
for bulk modulus $K_T$, and at 10000~K $\Delta$ is 29~\% and 3~\% for
same pressures.

\begin{table}[h]
\begin{tabular}{c r r r r r r r r}
%\hline
  T (K) & $c_s$ & $c_{44}$ & $c_{11}^T$ & $c_{11}^S$ &
$c_{12}^T$ & $c_{12}^S$ & $K_T$ & $K_S$ \\
\hline
   0 & 44.05 & 70.26 & 249.68 & 249.68 & 161.59 & 161.59 & 190.95 & 190.95 \\
 947 & 43.58 & 56.24 & 221.71 & 233.41 & 134.55 & 146.27 & 163.61 & 175.33 \\
2053 & 38.61 & 59.39 & 189.66 & 211.95 & 112.44 & 134.72 & 138.18 & 160.46 \\
3000 & 31.57 & 62.48 & 162.55 & 192.74 &  99.41 & 129.60 & 120.46 & 150.65 \\
3947 & 28.72 & 64.47 & 142.15 & 179.77 &  84.70 & 122.32 & 103.85 & 141.47 \\
5052 & 18.34 & 59.70 & 107.55 & 154.15 &  70.87 & 117.47 &  83.10 & 129.70 \\
6000 &  0.31 & 62.05 &  64.91 & 120.03 &  64.29 & 119.40 &  64.50 & 119.61 \\
%\hline
\end{tabular}
\vspace*{0.2cm}
\caption{\small The elastic constants for bcc Tantalum at various
temperatures. All elastic constants units are GPa.}
\label{table:econsther}
\end{table}

The shear moduli $c_s$ and $c_{44}$ and adiabatic bulk modulus $K_S$ agree
well with the ultrasonic experimental data~\cite{walker} up to 3000~K at
zero pressure (Fig.~\ref{fig:econsther}). We find that all three moduli are
primarily functions of volume, and thermal effects at constant volume are
quite small except at the highest pressures. There is some softening of $c_s$
with increasing temperature for all pressures. $c_{44}$, shows a slight
softening at the zero pressure with increasing temperature, but they are
rather flat for other pressures except than very high pressures.
The adiabatic bulk modulus $K_S$ also softens slightly with temperature at
low pressures but becomes flat with increasing pressure.

\begin{figure}
\centerline{\psfig{figure=anist.ps,angle=-90,width=80mm}}
%\vspace*{1cm}
\caption{\small The anisotropy ratio of elastic constants of Ta as a
function of temperature at different pressures from 0~GPa to 400~GPa with
50~GPa interval.}
\label{fig:anist}
\end{figure}

The anisotropy ratio, $A=c_{44}/c_s$, is presented as a function of
temperature for various pressures at Fig.~\ref{fig:anist}.
$A$ increases with increasing temperature at all pressures, but less
drastically at high pressures. At lower pressures, this increase is
divergent after certain temperature, since the softening of $c_s$ is
large enough and it approaches zero. The reversal of the slope of $A$
with pressure shifts to higher pressures with increasing temperature due
to thermal expansivity, and occurs at a fixed volume.

Sound velocities are related to the elastic constants by the Christoffel
equation~\cite{christoffel}
\begin{equation}
(c_{ijkl} n_j n_k - \rho v^2 \delta_{ij}) u_i = 0,
\end{equation}
where $c_{ijkl}$ is the elastic constants tensor, $\vec{n}$ is the
propagation direction, $\vec{u}$ is the polarization vector and
$v$ is the velocity. Our elastic constants are those appropriate for the
equations of motion under hydrostatic reference stress.\cite{barron}
For $[110]$ wave propagation direction in a cubic lattice, the
longitudinal mode is
\begin{equation}
\rho v^2 = (c_{11} + c_{12} + 2c_{44})/2
\end{equation}
and two transverse modes are
\begin{equation}
\rho v^2 = c_{44}
\end{equation}
and
\begin{equation}
\rho v^2 = (c_{11} - c_{12})/2 = c_s
\end{equation}
polarized along $[001]$ and $[1\bar{1}0]$ direction respectively.
For polycrystalline sample, the average isotropic shear modulus $G$ can
be determined from single crystal elastic constants according to the
Voigt-Reuss-Hill scheme~\cite{vrhaveg}, and the isotropically averaged
aggregate velocities are given by
\begin{eqnarray}
v_P & = & \left( (K+4/3G)/\rho\right)^{1/2} \\
v_S & = & (G/\rho)^{1/2} \\
v_B & = & (K/\rho)^{1/2}
\end{eqnarray}
where $v_P$, $v_S$, and $v_B$ are the compressional, shear and bulk
sound velocities. The sound velocities of Ta along the Hugoniot calculated
from elastic constants are shown in Fig.~\ref{fig:sound}, and are
compared with the shock sound velocity data from Brown {\it et al.}\cite{hug2}
As seen in Fig.~\ref{fig:sound}, there is excellent agreement with
shock data. The calculated compressional velocity $v_P$ agrees very well
with experimental data up to 200~GPa, and then after 300~GPa the bulk
velocity $v_B$ matches the data well. This is because the shocked
solid melts around 300~GPa, so the liquid velocity might be represented
by $V_B$. The deviation between 200~GPa and 300~GPa is probably due to
premelting effects.

\begin{figure}
\centerline{\psfig{figure=sound.ps,angle=-90,width=80mm}}
%\vspace*{1cm}
\caption{\small Sound velocities of Ta along the Hugoniot calculated 
from elastic constants. Solid lines are the longitudinal and two
transverse sound velocities in $[110]$ direction from single crystal
elastic constants. The polarization of the shear waves are given
in brackets. The isotropic aggregate velocities are shown by dashed
lines. $v_P$, $v_B$, and $v_S$ are the compressional, bulk, and shear
sound velocities.
Filled dots are the shock data from
Brown and Shaner~\protect\cite{hug2}.}
\label{fig:sound}
\end{figure}

In conclusion, the elasticity of bcc Ta is investigated from first
principles for pressures up to 400~GPa and temperatures up to 10000~K.
The calculated static elastic constants are in good agreement with
available ultrasonic and SAX experimental data. The shear
modulus $c_{44}$ and the anisotropy ratio $A$ change behaviour with
increasing pressure around 150~GPa. Although, the shear modulus $c_s$
softens with increasing temperature at all pressures, $c_{44}$ and
$K_S$ soften with temperature at low pressures but then they are rather
flat at higher pressures. The main effect of temperature for the
thermoelasticity of Ta is due to thermal expansivity. The calculated
sound velocities along the Hugoniot shows an excellent agreement with
shock-wave experimental data.

{\bf Acknowledgments}
We thank G. Steinle-Neumann, L. Stixrude, and E.  Wasserman for
helpful discussions. This work was supported by DOE ASCI/ASAP
subcontract B341492 to Caltech DOE W-7405-ENG-48.
Computations were performed on the Cray SV1 at the Geophysical Laboratory,
supported by NSF grant EAR-9975753 and the W.\ M.\ Keck Foundation.

\begin{references}

\bibitem{cynn} %Ta EOS experiment
H. Cynn and C.S, Yoo,
{\sl Phys. Rev. B} {\bf 59}, 8526 (1999).

\bibitem{hug1}  % Ta hugoniot
A.C. Mitchell and W.J. Nellis,
{\sl J. Appl. Phys.} {\bf 52}, 3363 (1981).

\bibitem{soderlind} % Livermore (LMTO) ta paper
P. S{\"o}derlind and J.A. Moriarty,
{\sl Phys. Rev. B} {\bf 57}, 10340 (1998).

\bibitem{tacohen} %Thermal equation of state of tantalum. 
R.E. Cohen, and O. G\"ulseren,
{\sl Phys. Rev. B} {\bf 63} 224101 (2001).

\bibitem{vinet} % Vinet universal EOS
P. Vinet, J. Ferrante, J.R. Smith and J.H. Rose,
{\sl J. Phys. C: Solid State Phys.} {\bf 19}, L467 (1986);
P. Vinet, J.H. Rose, J. Ferrante and J.R. Smith,
{\sl J. Phys.: Condens. Matter} {\bf 1}, 1941 (1987).

\bibitem{ron1}  % EOS paper
R.E. Cohen, O. G\"{u}lseren and R.J. Hemley,
{\sl Am. Mineral.} {\bf 85}, 338 (2000).

\bibitem{barron} %Initial stress correction
T.H.K. Barron, and M.L. Klein,
{\sl Proc. Phys. Soc.} {\bf 85} 523 (1965).

\bibitem{fecons}
G. Steinle-Neumann, L. Stixrude, R.E. Cohen
and O. G\"{u}lseren,
{\sl Nature} {\bf 413}, 57 (2001).

\bibitem{wass}
E. Wasserman, L. Stixrude and R.E. Cohen,
{\sl Phys. Rev. B} {\bf 53}, 8296 (1996).

\bibitem{iron} %hcp iron
D.Alfe, G.D. Price, and M.J. Gillan,
{\sl Phys. Rev. B} {\bf 64} 045123 (2001).

\bibitem{temdft} % finite temperature DFT
N.D. Mermin,
{\sl Phys. Rev.} {\bf 137}, A1441 (1965).

\bibitem{weikrakauer}
S.H. Wei and H. Krakauer,
{\sl Phys. Rev. Lett.} {\bf 55},  1200-1203  (1985).

\bibitem{singh}  % D. Singh's LAPW book
D.J. Singh,
{\sl Planewaves, Pseudopotentials, and the LAPW Method}
(Kluwer Academic Publishers, Boston, 1994).

\bibitem{gga-pbe} % GGA-PBE
J.P. Perdew, K. Burke and M. Ernzerhof,
{\sl Phys. Rev. Lett.} {\bf 77}, 3865 (1996).

\bibitem{monpack} % Monk-Pack special k-points
H.J. Monkhorst and J.D. Pack,
{\sl Phys. Rev. B} {\bf 13}, 5188 (1976).

\bibitem{cellmodel} % Cell model papers
A.C. Holt, W.G. Hoover, S.G. Gray and D.R. Shortle,
{\sl Physica} {\bf 49}, 61 (1970);
F.H. Ree and A.C. Holt,
{\sl Phys. Rev. B} {\bf 8}, 826 (1973);
K. Westera and E.R. Cowley,
{\sl Phys. Rev. B} {\bf 11}, 4008 (1975);
E.R. Cowley, J. Gross, Z. Gong and G.K. Horton,
{\sl Phys. Rev. B} {\bf 42}, 3135 (1990).

\bibitem{mix-oguz}  % Implementation of mixed basis method
O. G\"{u}lseren, D.M. Bird and S.E. Humphreys,
{\sl Surf. Sci.} {\bf 402-404}, 827 (1998).

\bibitem{grigoriev} 
{\it Handbook of Physical Quantities}, Ed. I.S.
Grigoriev and E.Z. Meilikhov, CRC Press, NY (1997).

\bibitem{sonali} % vacancy formation energy
S. Mukherjee, R.E. Cohen and  O. G\"{u}lseren,
{\sl submitted to Phys. Rev. B} 

\bibitem{ldazun} % LDA parametrization
J.P. Perdew and A. Zunger,
{\sl Phys. Rev. B} {\bf 23}, 5048 (1981).

\bibitem{tromartins} %Cu pseudopotential
N. Troullier and J.L. Martins,
{\sl Phys. Rev. B} {\bf 43}, 1993 (1991).

\bibitem{nonlinearcore} % non-linear core correction
S.G. Louie, S. Froyen and M.L. Cohen,
{\sl Phys. Rev. B} {\bf 26}, 1738 (1982).

\bibitem{spdir} %integration with special directions
A. Bansil,
{\sl Solid State Commun.} {\bf 16}, 885 (1975);
R. Prasad and A. Bansil,
{\sl Phys. Rev. B} {\bf 21}, 496 (1980);
W.R. Fehlner, S.B. Nickerson and S.H. Vosko,
{\sl Solid State Commun.} {\bf 19}, 83 (1976).

\bibitem{katahara1}
K.W. Katahara, M.H. Manghnani and E.S. Fisher,
{\sl J. Phys. F: Metal Phys.} {\bf 9}, 773 (1979).

\bibitem{katahara2}
K.W. Katahara, M.H. Manghnani and E.S. Fisher,
{\sl J. Appl. Phys.} {\bf 47}, 434 (1976).

\bibitem{cynnecons} %Ta Econs experiment
H. Cynn and C.S, Yoo,
{\sl Submitted.}; H. Cynn and C.S, Yoo,
{\sl Bull. Am. Phys. Soc.} {\bf 44}, 1288 (1999).

\bibitem{davies} %Adiabatic elastic constants
G.F. Davies,
{\sl J. Phys. Chem. Solids} {\bf 35}, 1513 (1974).

\bibitem{walker} %Ta high temperature econs experiment
E. Walker and P. Bujard,
{\sl Solid State Comm.} {\bf 34}, 691 (1980).

\bibitem{christoffel} %Christoffel equation
E.B. Christoffel,
{\sl Ann. Mat. Pura Appl.} {\bf 8}, 193 (1877).

\bibitem{vrhaveg} %Voigt-Reuss-Hill average
R.W. Hill,
{\sl Proc. Phys. Soc. London} {\bf 65 A}, 349 (1952).

\bibitem{hug2}  % Ta hugoniot, sound velocities
J.M. Brown and J.W. Shaner,
{\it Shock Waves in Condensed Matter-1983},
edited by J.R. Asay, R.A. Graham and G.K. Straub,
pp 91-94, Elsevier, New York (1984).

\end{references}

%\clearpage

\end{document}

cond-mat0111227.extracted_bib
\bibitem{cynn} %Ta EOS experiment
H. Cynn and C.S, Yoo,
{\sl Phys. Rev. B} {\bf 59}, 8526 (1999).


\bibitem{hug1}  % Ta hugoniot
A.C. Mitchell and W.J. Nellis,
{\sl J. Appl. Phys.} {\bf 52}, 3363 (1981).


\bibitem{soderlind} % Livermore (LMTO) ta paper
P. S{\"o}derlind and J.A. Moriarty,
{\sl Phys. Rev. B} {\bf 57}, 10340 (1998).


\bibitem{tacohen} %Thermal equation of state of tantalum. 
R.E. Cohen, and O. G\"ulseren,
{\sl Phys. Rev. B} {\bf 63} 224101 (2001).


\bibitem{vinet} % Vinet universal EOS
P. Vinet, J. Ferrante, J.R. Smith and J.H. Rose,
{\sl J. Phys. C: Solid State Phys.} {\bf 19}, L467 (1986);
P. Vinet, J.H. Rose, J. Ferrante and J.R. Smith,
{\sl J. Phys.: Condens. Matter} {\bf 1}, 1941 (1987).


\bibitem{ron1}  % EOS paper
R.E. Cohen, O. G\"{u}lseren and R.J. Hemley,
{\sl Am. Mineral.} {\bf 85}, 338 (2000).


\bibitem{barron} %Initial stress correction
T.H.K. Barron, and M.L. Klein,
{\sl Proc. Phys. Soc.} {\bf 85} 523 (1965).


\bibitem{fecons}
G. Steinle-Neumann, L. Stixrude, R.E. Cohen
and O. G\"{u}lseren,
{\sl Nature} {\bf 413}, 57 (2001).


\bibitem{wass}
E. Wasserman, L. Stixrude and R.E. Cohen,
{\sl Phys. Rev. B} {\bf 53}, 8296 (1996).


\bibitem{iron} %hcp iron
D.Alfe, G.D. Price, and M.J. Gillan,
{\sl Phys. Rev. B} {\bf 64} 045123 (2001).


\bibitem{temdft} % finite temperature DFT
N.D. Mermin,
{\sl Phys. Rev.} {\bf 137}, A1441 (1965).


\bibitem{weikrakauer}
S.H. Wei and H. Krakauer,
{\sl Phys. Rev. Lett.} {\bf 55},  1200-1203  (1985).


\bibitem{singh}  % D. Singh's LAPW book
D.J. Singh,
{\sl Planewaves, Pseudopotentials, and the LAPW Method}
(Kluwer Academic Publishers, Boston, 1994).


\bibitem{gga-pbe} % GGA-PBE
J.P. Perdew, K. Burke and M. Ernzerhof,
{\sl Phys. Rev. Lett.} {\bf 77}, 3865 (1996).


\bibitem{monpack} % Monk-Pack special k-points
H.J. Monkhorst and J.D. Pack,
{\sl Phys. Rev. B} {\bf 13}, 5188 (1976).


\bibitem{cellmodel} % Cell model papers
A.C. Holt, W.G. Hoover, S.G. Gray and D.R. Shortle,
{\sl Physica} {\bf 49}, 61 (1970);
F.H. Ree and A.C. Holt,
{\sl Phys. Rev. B} {\bf 8}, 826 (1973);
K. Westera and E.R. Cowley,
{\sl Phys. Rev. B} {\bf 11}, 4008 (1975);
E.R. Cowley, J. Gross, Z. Gong and G.K. Horton,
{\sl Phys. Rev. B} {\bf 42}, 3135 (1990).


\bibitem{mix-oguz}  % Implementation of mixed basis method
O. G\"{u}lseren, D.M. Bird and S.E. Humphreys,
{\sl Surf. Sci.} {\bf 402-404}, 827 (1998).


\bibitem{grigoriev} 
{\it Handbook of Physical Quantities}, Ed. I.S.
Grigoriev and E.Z. Meilikhov, CRC Press, NY (1997).


\bibitem{sonali} % vacancy formation energy
S. Mukherjee, R.E. Cohen and  O. G\"{u}lseren,
{\sl submitted to Phys. Rev. B} 


\bibitem{ldazun} % LDA parametrization
J.P. Perdew and A. Zunger,
{\sl Phys. Rev. B} {\bf 23}, 5048 (1981).


\bibitem{tromartins} %Cu pseudopotential
N. Troullier and J.L. Martins,
{\sl Phys. Rev. B} {\bf 43}, 1993 (1991).


\bibitem{nonlinearcore} % non-linear core correction
S.G. Louie, S. Froyen and M.L. Cohen,
{\sl Phys. Rev. B} {\bf 26}, 1738 (1982).


\bibitem{spdir} %integration with special directions
A. Bansil,
{\sl Solid State Commun.} {\bf 16}, 885 (1975);
R. Prasad and A. Bansil,
{\sl Phys. Rev. B} {\bf 21}, 496 (1980);
W.R. Fehlner, S.B. Nickerson and S.H. Vosko,
{\sl Solid State Commun.} {\bf 19}, 83 (1976).


\bibitem{katahara1}
K.W. Katahara, M.H. Manghnani and E.S. Fisher,
{\sl J. Phys. F: Metal Phys.} {\bf 9}, 773 (1979).


\bibitem{katahara2}
K.W. Katahara, M.H. Manghnani and E.S. Fisher,
{\sl J. Appl. Phys.} {\bf 47}, 434 (1976).


\bibitem{cynnecons} %Ta Econs experiment
H. Cynn and C.S, Yoo,
{\sl Submitted.}; H. Cynn and C.S, Yoo,
{\sl Bull. Am. Phys. Soc.} {\bf 44}, 1288 (1999).


\bibitem{davies} %Adiabatic elastic constants
G.F. Davies,
{\sl J. Phys. Chem. Solids} {\bf 35}, 1513 (1974).


\bibitem{walker} %Ta high temperature econs experiment
E. Walker and P. Bujard,
{\sl Solid State Comm.} {\bf 34}, 691 (1980).


\bibitem{christoffel} %Christoffel equation
E.B. Christoffel,
{\sl Ann. Mat. Pura Appl.} {\bf 8}, 193 (1877).


\bibitem{vrhaveg} %Voigt-Reuss-Hill average
R.W. Hill,
{\sl Proc. Phys. Soc. London} {\bf 65 A}, 349 (1952).


\bibitem{hug2}  % Ta hugoniot, sound velocities
J.M. Brown and J.W. Shaner,
{\it Shock Waves in Condensed Matter-1983},
edited by J.R. Asay, R.A. Graham and G.K. Straub,
pp 91-94, Elsevier, New York (1984).

