[gpumd] Re: Can I extract specific heat information from dos.out file?

  • From: Bruce Fan <brucenju@xxxxxxxxx>
  • To: gpumd@xxxxxxxxxxxxx
  • Date: Thu, 21 Nov 2019 18:04:03 +0200

Hi,

Yes you can. This is a piece of cake.

Check the updated Matlab script here (
https://github.com/brucefan1983/GPUMD/blob/master/examples/gpumd/ex2/plot_results.m)
which can reproduce the figure below:

[image: untitled.png]

On Mon, Nov 18, 2019 at 3:01 PM "소순성" <soonsung2001@xxxxxxxxxx> wrote:

Dear Zheyoung Fan,


Thank you  Zheyong Fan for using this good GPUMD program.

Today, I have a question.


I wonder if I can extract the specific heat information from dos.out file?


because I'm poor at coding, I fixed plot.m files in "ex#" folders in gpumd
package so that it contains:


clear; close all; font_size=10;

load mvac.out;

N=250; % number of correlation steps

M_v=length(mvac)/N % number of independent functions


load dos.out;

nu=dos(1:N,1)/2/pi; % nu = omega/2/pi

M_d=length(dos)/N

dos_x=mean(reshape(dos(:,2),N,M_d),2);

dos_y=mean(reshape(dos(:,3),N,M_d),2);

dos_z=mean(reshape(dos(:,4),N,M_d),2);


% calculated parameters

dt=4;

dt_in_ps = dt/1000;       % ps

time_in_ps=(0:N-1)*dt_in_ps;


atoms=15228;

Na=6.0221409e+23;

kb=1.3806488e-23;        %(J/K)

kbev=8.617343e-5;        %(eV/K)

h=6.62607004081e-34;  %(J*s)

hbar=h/pi/2;

hev=4.1356678064e-15;  %(eV*s)

T=300;

const_nu=h/(kb)/2/pi/T;

const_nu=hbar/kb/T;

const_nuev=hev/(kbev)/2/pi/T;



%Calculation of Cv according to axis


for n=1:length(nu)

 cx(n)=
2*pi*nu(n)*dos_x(n)*((const_nu)*2*pi*nu(n)).^2.*exp((const_nu)*2*pi*nu(n))./(exp((const_nu)*2*pi*nu(n))-1).^2.;

 cy(n)=
2*pi*nu(n)*dos_y(n)*((const_nu)*2*pi*nu(n)).^2.*exp((const_nu)*2*pi*nu(n))./(exp((const_nu)*2*pi*nu(n))-1).^2.;

 cz(n)=
2*pi*nu(n)*dos_z(n)*((const_nu)*2*pi*nu(n)).^2.*exp((const_nu)*2*pi*nu(n))./(exp((const_nu)*2*pi*nu(n))-1).^2.;

end



Cvx=(2*pi)*cx*nu;

Cvy=(2*pi)*cy*nu;

Cvz=(2*pi)*cz*nu;


nureshape=reshape(dos(1:N,1),1,N);

normalization_of_specific_heat1_x=trapz(nureshape,cx)

normalization_of_specific_heat1_y=trapz(nureshape,cy)

normalization_of_specific_heat1_z=trapz(nureshape,cz)



% An important check is that PDOS should be normalized to about 1/2

normalization_of_pdos_x=trapz(nu,dos_x)

normalization_of_pdos_y=trapz(nu,dos_y)

normalization_of_pdos_z=trapz(nu,dos_z)


But I think it is uncorrect and I tried many times but It was failed.

If it is possible, how can I do?

How can I get specific heat value from dos.out?











PNG image

Other related posts: