24小时热门版块排行榜

返回列表

雁儿霏霏

木虫 (正式写手)

应助: 2 (幼儿园)
金币: 2645.4
帖子: 883
在线: 160.1小时
虫号: 934537

[交流] 【求助】紧急求教读取CHGCAR的matlab脚本，谢谢

现在需要深入研究vasp的CHGCAR，希望用matlab读取它，但对矩阵的读取一直没查到相关资料，求教哪位有读能提供取CHGCAR的matlab小脚本，谢谢。

回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

继续求助Matlab绘图！已经有8人回复
弱弱问一下：stereochemical activity 是什么意思？谢谢已经有5人回复
紧急求助翻译一篇德文文献~翻译其中合成的一部分~急求~谢谢~ 已经有4人回复
matlab 计算过程中工作空间变量保存及读取求解决内存不足的方法已经有5人回复
帮翻译几段文献，请不要在线翻译，谢谢已经有2人回复
用matlab统计图片中某种颜色成分占总图片中所有颜色成分的比率已经有4人回复
【求助】弱弱的问一下Matlab数据提取的问题已经有5人回复
【求助】用matlab编程，需要将数据自检分类（高手进来指点下）已经有5人回复
【求助】CHGCAR用什么来读取已经有10人回复
【求助】matlab如何只读取图像的一部分？已经有11人回复
【求助】有哪位大侠做过次氯酸叔丁酯啊？帮帮我啊！谢谢已经有6人回复
【求助】matlab读取fortran输出的txt文件【已解决】已经有3人回复
【求助】求教matlab解非线性方程组已经有9人回复
【求助】Weickert的各向异性扩散方程的滤波方法的matlab程序代码已经有8人回复

» 抢金币啦！回帖就可以得到:

查看全部散金贴

1楼 2011-01-15 10:33:48

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

liangab_234620

木虫 (小有名气)

应助: 0 (幼儿园)
金币: 1341.2
帖子: 131
在线: 145.5小时
虫号: 461188

★ ★
sunyang1988(金币+2): 谢谢帮助。这样直接粘贴可能有格式错误，如果方便请上传到网盘 2011-01-16 14:10:22
雁儿霏霏(金币+20): 谢谢大虾，能给提供个程序文件吗？粘贴有一些格式错误。邮箱：shangyan2009@gmail.com ，谢谢！ 2011-01-16 19:37:24

% Read in CHG file to 3D matrix.

fid = fopen('CHG','r');

% Read POSCAR part of CHG file:
sysname = fgetl(fid);
lconst = fscanf(fid,'%f',1);
basisvec = fscanf(fid,'%f',[3,3])';
temp1 = fgetl(fid);
temp2 = fgetl(fid);
species = sscanf(temp2,'%i');
Natoms = sum(species);
temp3 = fgetl(fid);
atompos_f3 = fscanf(fid,'%f',[3,Natoms])';

volume = abs(dot(basisvec(1,

,cross(basisvec(2,

,basisvec(3,

)))*lconst^3;

% Read charge density data. (i,j,k) value in i+(j+k*Ny)*Nx, with coordinate
% system defined by basis vectors!:
gridsize = fscanf(fid,'%i',3);
Nx = gridsize(1);
Ny = gridsize(2);
Nz = gridsize(3);

temp4 = fscanf(fid,'%f',[Nx*Ny*Nz,1]);

density_sub_h1_surf=reshape(temp4,Nx,Ny,Nz)/volume;

% % Magnitization (if present):
% temp5 = fscanf(fid,'%i',3);
% temp6 = fscanf(fid,'%f',[Nx*Ny*Nz,1]);
% magn_lowacc_t2_N�蕆eshape(temp6,Nx,Ny,Nz)/volume;

% minval = min(min(min(density)));
% maxval = max(max(max(density)));

% For rectangular unit cell, unit axes:
% Plot isosurface:
% prect = patch(isosurface(density,7.5), 'FaceColor', 'red', 'EdgeColor', 'none');
% isonormals(density,prect);
% %axis tight;
% camlight; lighting phong;

% Non-rectangular (general) unit cell:
% Need to created position 3D matrices
xtemp = zeros(Nx*Ny*Nz,1);
ytemp = zeros(Nx*Ny*Nz,1);
ztemp = zeros(Nx*Ny*Nz,1);

for k = 0:Nz-1
for j = 0:Ny-1
      for i = 0:Nx-1
         xtemp(1+i+(j+k*Ny)*Nx) = i/Nx*basisvec(1,1)+j/Ny*basisvec(2,1)+k/Nz*basisvec(3,1);
         ytemp(1+i+(j+k*Ny)*Nx) = i/Nx*basisvec(1,2)+j/Ny*basisvec(2,2)+k/Nz*basisvec(3,2);
         ztemp(1+i+(j+k*Ny)*Nx) = i/Nx*basisvec(1,3)+j/Ny*basisvec(2,3)+k/Nz*basisvec(3,3);
      end
end
end

X = reshape(xtemp,Nx,Ny,Nz)*lconst;
Y = reshape(ytemp,Nx,Ny,Nz)*lconst;
Z = reshape(ztemp,Nx,Ny,Nz)*lconst;

z = squeeze(Z(1,1,

);
ave_z_h1 = squeeze(sum(sum(density_sub_h1_surf,1),2))/(Nx*Ny); %(sqrt(2)*Nx*Ny)*lconst^2;

figure
%subplot(211)
plot(z,ave_z_h1,'k')

% figure
% p = patch(isosurface(X,Y,Z,density,7.5), 'FaceColor', 'red', 'EdgeColor', 'none');
% isonormals(density,p);
% % maxx = max(basisvec(:,1))*lconst;
% % minx = min(basisvec(:,1))*lconst;
% % maxy = max(basisvec(:,2))*lconst;
% % miny = min(basisvec(:,2))*lconst;
% % maxz = max(basisvec(:,3))*lconst;
% % minz = min(basisvec(:,3))*lconst;
% % axis([minx maxx miny maxy minz maxz])
% camlight; lighting phong;

% Extract charge density along certain lines:

------------------------------