前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >PVOX-自定义函数readoutput分析

PVOX-自定义函数readoutput分析

作者头像
巴山学长
发布2020-05-08 11:11:10
4340
发布2020-05-08 11:11:10
举报
文章被收录于专栏:巴山学长巴山学长

“庖丁为文惠君解牛,手之所触,肩之所倚,足之所履,膝之所踦,砉然向然,奏刀騞然,莫不中音。合于《桑林》之舞,乃中《经首》之会”。

牛虽难解,十载不费一刀,徒心中有牛。复杂工程不要慌,看过冷水带你一步一步解析PVOX工具包,本期先看看自定义函数readoutput函数的构建。

代码语言:javascript
复制
[data, wfn, D_valid, W_valid] = readOutput(outName, datName)

(1)自定义函数readOutput()的作用是读取后缀为*.out和*.dat的两个文件,输出data、wfn、D_valid、W_valid对象;

(2)Data:的作用是提取*.out中的一些变量和对应的数据,重新储存在data对象中。

(3)wfn:是提取后缀为*.dat的文件中的数据,储存在wfn中。表征波函数

(4)D_valid、W_valid这两个量是用于监控Data、wfn过程环节是否出错而设置的。假设wfn计算没有问题则W_valid赋值为1,否则为0;

代码语言:javascript
复制
outName='parsec_grid0_4.out'
datName='parsec_grid0_4.dat'
D_valid = false;
W_valid = [0,0,0,0,0];
data = [];
wfn = [];

该步骤为导入两个文件。定义输出对象。需要注意的是: D_valid为单值,W_valid为多值,这是因为存储wfn数据过程中有多个子环节需要判断。

接下是try-catch-end结构

代码语言:matlab
复制
try   
    Statements1
catch    
    Statements2
end

该结构的含义是尝试执行命令。先执行try层下的语句命令sta1,若是正常执行,则该结构功能结束,若是try层语句命令不能正常执行,则执行catch 层下的语句命令sta2,如果sta1、sta2皆不能正常执行,则跳过该结构代码,执行后续命令,不报错。这样来配合变量监控就可以从整体来把控代码有几处问题了,而不是代码以有错,就停止判断,无法了解后续代码情况。

代码语言:matlab
复制
h = waitbar(0,'Opening parsec.out');
[str,maxsize,endian] = computer;
fName = fopen('pTemp/fname.dat','wt');
fName = fopen('pTemp/fname.dat','wt');
fwrite(fName,datName);
fclose(fName);
waitbar(0.1,h);

该部分的语句的含义是新建一个文件,将后缀为*.out和*.dat的两个文件的路径写入该文件中,以便后续使用。

原始代码的第一个错误就在这里出现了。我们以截图的形式比较两段代码的差别。

过冷水的修改代码中多了两行命令,

代码语言:matlab
复制
fprintf(fName,'\n');
fwrite(fName,datName);

What you want do!过冷水总会为自己擅自修改代码编造出合理的理由,我们结合这两段代码来看另一段代码。

如果不对matlab中fname.dat写入文件时进行适当修改,下图的代码2就会报错。代码1、2是python码,不了解python的加入我们python课程,懂一点python让我们的距离更加近一点!在此说明一下:

代码语言:matlab
复制
fname = file('./pTemp/fname.dat', 'r')
fileName = fname.readline()

目的为打开fname.dat文件,并读取第一行的内容。我们们将*.out文件路径写入fname.dat中,让其读取,没毛病!

代码语言:matlab
复制
fname = file('./pTemp/fname.dat', 'r')
fileName = fname.readline()
fileName = fname.readline()
switchEndian = fname.readline()
inFile = file(fileName, 'rb')

目的为打开fname.dat文件,并读取第二行、第三行的内容。什么时候往fname.dat文件中第二行和第三行写入数据了?能不报错吗?这就是为什么添加写入内容的原因。这里改写需要注意一下各种error!

想要给读者展示错误代码报错,容易的很,在此过冷水只展示了几种没有明显错误,但实际是错的案例。过冷水只是想往文本里多添一条绝对路径,需要注意的点就有这么多,可见代码的编写细节很多,不断学习才能够完善编程知识。需要你精通matlab的跟着过冷水一行一行看代码!

代码语言:matlab
复制
if strncmp(str,'PCWIN',5)        
    try   % Parsec Abridged Output            
        [s,w] = dos('del /Q pTemp\main.dat');            
        !"./pTools/PARSEC_PARSER_ABRIDGED.py"        
    catch % Parsec Un-Abridged Output            
        !"./pTools/PARSEC_PARSER.py"        
    end  
else:   
    [s,w] = unix('rm -f pTemp/main.dat');

新知识[!"fName"]该格式是在matlab中打开文件的命令。

该段代码是判断我们运行matlab的系统是什么,不同的系统用不同的方式打开*.py文件。继续接着走:

代码语言:javascript
复制
inMain = load('pTemp/main.dat');   
    data.timeDiag = inMain(:,1);    
    data.FermiLevel = inMain(:,2);    
    data.timeHart = inMain(:,3);    
    data.ChargeDmin = inMain(:,4);    
    data.ChargeDmax = inMain(:,5);    
    data.Eig_Energy = inMain(:,6);    
    data.H_Energy = inMain(:,7);    
    data.Int_Vxc_rho = inMain(:,8);    
    data.Exc = inMain(:,9);    
    data.E_I_Energy = inMain(:,10);    
    data.I_I_Energy = inMain(:,11);    
    data.Enew_Eold = inMain(:,12);     
    data.Tot_Energy = inMain(:,13);    
    data.Energy_atom = inMain(:,14);    
    data.SREpot = inMain(:,15);
    data.ChrgWghtPot = inMain(:,16);

这段代码的目的为从main.dat(运行*.py文件得到)中导入数据,储存到data对象中。这里再次get一个新知识点。以前过冷水就只知道:

代码语言:matlab
复制
data.timeDiag = inMain(:,1)
data= inMain(:,1)

现在变了一种赋值形式:

代码语言:matlab
复制
>>  inMain = load('pTemp/main.dat');    
    data.timeDiag = inMain(:,1)
data =     
    timeDiag: [9x1 double]

Dat.char这种赋值方式很有用的,可以让我们很容易储存不同类型的数据到一个对象中,直接表示数据类型,不需要做其它操作,类似于元包组。

代码语言:matlab
复制
fOptions = fopen('pTemp/misc.dat');
opt = fgetl(fOptions);
data.num_states = str2num(strtok(opt));    
    opt = fgetl(fOptions);    
    data.num_atoms = str2num(strtok(opt));    
    iter = size(data.timeDiag);    
    data.num_iterations = iter(1);

这段代码有三个新函数:fgetl()、str2num();

fgetl():该函数的目的是每次运行从misc.dat文件中获取一行内容。fgetl()有记忆功能。这里要注意运行次数,稍微运行次数出错会导致写入的数据对不上号,其实这里如果可以用正则匹配或者关键字定位行就不容易出错了;

strtok():函数的含义是从字符串中找出数值字符串;

str2num:函数的目的是将字符串转化为数值。

代码语言:matlab
复制
inEigs = load('pTemp/eigenvalues.dat');    
    m = data.num_states;    
    n = size(inEigs(:,1));    
    n = n(1);    
    n = n/m;    
    data.State = reshape(inEigs(:,1),m,n);    
    data.EigenvalueRY = reshape(inEigs(:,2),m,n);    
    data.EigenvalueEV = reshape(inEigs(:,3),m,n);    
    data.Occup = reshape(inEigs(:,4),m,n);    
    data.Repr = reshape(inEigs(:,5),m,n)

这段代码就没有什么新意了,打开一个文件读取文件内容,并用一定规则重构元素排列。Reshape()、size()都是常用命令。

代码语言:matlab
复制
inForces = load('pTemp/forces.dat');    
    forceIon = reshape(inForces(1,:),3,data.num_atoms);    
    data.forceIon = transpose(forceIon);    
    forceSelfC = reshape(inForces(2,:),3,data.num_atoms);    
    data.forceSelfC = transpose(forceSelfC);    
    forceCoreC = reshape(inForces(3,:),3,data.num_atoms);    
    data.forceCoreC = transpose(forceCoreC);

代码调试不通,困扰过冷水的第二处错误就此出现。load()加载文件、reshape重构数据结构、transpose()函数没见过?调用出错?

函数transpose():看上去定义又长,好像很重要的样子,过冷水一查:

B=transpose(A)=A’结束!

既然语句看上去没有问题,那么会有什么问题呢?这里过冷水怎么都没有想到这个坑,逻辑错误!什么意思呢?你要加载这个文件,首先的有这个文件吧!没有文件怎么加载。如果没有这个文件那就是error呗!

解释的有点费劲。*.out生成文件一般默认不包含forceIon、 forceSelfC、forceCoreC数据。然而这里的程序包默认*.out含有相应数据,不做有无的判断!这是程序设计的漏洞,因为这个原因,*.py程序编写的也有问题。过冷水的解决办法是用NaN 填充相应数据。这里可以让程序运行正确,我们不考虑科学合理性的问题。

代码语言:matlab
复制
inForces2 = load('pTemp/forces2.dat');    
    forceLocal = reshape(inForces2(1,:),6,data.num_atoms);    
    data.forceLocal = transpose(forceLocal);    
    forceNonLcl = reshape(inForces2(2,:),6,data.num_atoms);    
    data.forceNonLcl = transpose(forceNonLcl);    
    forceTotal = reshape(inForces2(3,:),6,data.num_atoms);    
    data.forceTotal = transpose(forceTotal);    
    forceTotal2 = reshape(inForces2(4,:),6,data.num_atoms);    
    data.forceTotal2 = transpose(forceTotal2);
    inForces3 = load('pTemp/forces3.dat');    
    data.dipole = reshape(inForces3(1,:),1,3);    
    data.netForce = reshape(inForces3(2,:),1,3);    
    data.netTorque = reshape(inForces3(3,:),1,3);
D_valid = true;

data数据写入完毕!该段代码没有什么好讲的,和前面的功能相似,注意D_valid = true命令!这里就相当于在做是否有完整执行data数据写入的代码。该段代码让过冷水的涨见识了,简单的写数据的操作就几百行代码,可唬人了。经过冷水一番解刨,不过如此!

接下来回是wfn数据的写入,这里有涉及一个新自定义函数

代码语言:matlab
复制
[wfn, pot, rho, wavedata, XYZ] = wfnread(datName, h)

这真他*的是一个牵一发,动全身的工程。该完整函数程序比较复杂,初次学习涉及到的新函数比较多,需讲解内容较多,本期只讲data数据的提取。欲学后情,且下次再讲。

推荐指数:★★★★ (8/10分)

好不好用只有用了才知道!若觉得好,别忘分享给和您一样爱学习研究的小伙伴哦!

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-05-02,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 巴山学长 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档