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?