用MS的图形界面操作MOPAC计算的懒人包

Shannon

如题，利用material studio强大的建模能力快速 输入几何结构到mopac中， mopac计算完得到的几何结构重新输入到material studio里面。 结果立等可取，找一个过渡态只需要2秒钟。
例如输入结构：


输出结构


如果在输入结构中测量了键长键角，输出的结构中也会直接显示出键长键角，连点鼠标看键长的操作也可以省掉。
代码如下（MATLAB代码）：
主要执行文件

1. [atomsfetched,xsdfile]= xml2atoms('TSoutput.xsd');
   
2. [xyz_num,species] = atoms2xyz7name (atomsfetched);
   
3. filename='rhotbTS'
   
4. xyz2mopac_inpTS (filename,xyz_num,species);
   
5. dos(['mopac2012 ',filename,'.mop'])
   
6. %%
   
7. [species,xyz_out] = findgeommopac([filename,'.out']);
   
8. %output section
   
9. xyz_num2atoms (xyz_out,atomsfetched);% no need to output anything, atomsfetched is merely a reference.
   
10. xmlwrite( 'TSoutput.xsd',xsdfile) %将对象存入电脑

复制代码

用来读取从xsd文件中提取 需要的atom对象的程序：xsd读取代码用的是java的语法…读起来感觉可能会很奇怪。

1. function [atomsfetched,xsdfile]= xml2atoms(pathofxmkl)
   
2. xsdfile= xmlread(pathofxmkl);
   
3. root=xsdfile.getDocumentElement;
   
4. root.setAttribute('Version','6.0');
   
5. root.setAttribute('WrittenBy','Materials Studio 6.0');
   
6. tmpnode=root.getFirstChild;
   
7. nodeofinterest=tmpnode.getNextSibling;
   
8. 
   
9. for ix=0:nodeofinterest.getLength()-1% Java 是从0开始计数的！！！
   
10.     currentitem=nodeofinterest.item(ix);%
    
11.     if strcmp(char(currentitem.getClass),'class org.apache.xerces.dom.DeferredElementImpl')
    
12. if currentitem.getTagName=='Molecule'
    
13.     moleculeinfo=nodeofinterest.item(ix);
    
14. end
    
15.     end
    
16. end
    
17. 
    
18. ilist=1;
    
19. for iy=0:moleculeinfo.getLength()-1
    
20. 
    
21.       currentatomitem=moleculeinfo.item(iy);
    
22. 
    
23.         if strcmp(char(currentatomitem.getClass),'class org.apache.xerces.dom.DeferredElementImpl')
    
24.              if currentatomitem.getTagName=='Atom3d'
    
25.                
    
26. atomsfetched(ilist)=currentatomitem;
    
27. ilist=ilist+1;
    
28.              end
    
29.         end
    
30.   
    
31. end
    
32. end

复制代码

从atom对象中提取坐标信息：

1. function [xyz_num,species] = atoms2xyz7name (atomsfetched)
   
2. xyz_num=zeros(length(atomsfetched),3);
   
3. for iz=1:length(atomsfetched)% 改用matlab语法
   
4.    currentatom=atomsfetched(iz);
   
5.    xyz_java_str(iz)=currentatom.getAttribute('XYZ');
   
6.    xyz_matstr=char(xyz_java_str(iz)); % 转换java变量到matalb变量
   
7.    xyz_matcell=strsplit(xyz_matstr,',');%从逗号那儿分割整个字符串
   
8.    xyz_num(iz,:)= str2double(xyz_matcell);%转换字符形式为数字
   
9.      species(iz)=currentatom.getAttribute('Components');
   
10. end
    
11. species=char(species);
    
12. end

复制代码

将坐标信息写成mopac输入文件

1. function output = xyz2mopac_inp (filename,xyz_num,species)
   
2. mopacconfigure=['PM7 eps=78.54 TS  CHARGE=0 \n my mopac pre optimization \n Name card requried \n '];
   
3. 
   
4. xyzinput='';
   
5. [nspec,~]=size(xyz_num);
   
6. for ispec=1:nspec
   
7.     xyzinput = [xyzinput,species(ispec),' ',num2str(xyz_num(ispec,:),'%010.9f '),'\n'];
   
8. end
   
9. xyzinput=sprintf(xyzinput);
   
10. 
    
11. inputfilesetting_mopac=[mopacconfigure,xyzinput];
    
12. inputfile_mopac=sprintf(inputfilesetting_mopac);
    
13. 
    
14. fid1=fopen([filename,'.mop'],'w');
    
15. fprintf(fid1,'%s',inputfile_mopac)
    
16. fclose(fid1);
    
17. output=0;
    
18. end

复制代码

从mopac程序输出文件中提取坐标信息

1.   function [species,xyz] = findgeommopac(filename)
   
2. fid=fopen(filename,'r');
   
3. 
   
4. ix=1;
   
5. while 1
   
6.    iline= fgetl(fid);
   
7.    if iline==-1
   
8.     break
   
9.    end
   
10.     geominfo=regexp(iline,'CARTESIAN COORDINATES') ;
    
11.    if ~isempty(geominfo)
    
12.       filepositions(ix)=ftell(fid);
    
13.       ix=ix+1;
    
14.    end
    
15.    
    
16. end
    
17. outputgeompst=filepositions(end);
    
18. fseek(fid,outputgeompst,'bof'); % 鼠标位置转移到CARTESIAN COORDINATES的末尾，刚读取玩CARTESIAN COORDINATES这一行的时候
    
19. fgetl(fid);
    
20. iy=1;
    
21. while 1
    
22. currentdataline=  fgetl(fid) ;
    
23.    if isempty(currentdataline)
    
24.       break
    
25.   end
    
26.    data=strsplit(currentdataline);
    
27.    species(iy)=data{3};
    
28.    xyz(iy,1)=str2double(data{4});
    
29.    xyz(iy,2)=str2double(data{5});
    
30.    xyz(iy,3)=str2double(data{6});
    
31.    iy=iy+1;
    
32.    
    
33. end
    
34. 
    
35. fclose(fid);
    
36. end

复制代码

将读取到的坐标重新输入material studio

1. function useless =xyz_num2atoms (xyz_num,atomsfetched)
   
2. for iz=1:length(atomsfetched)   
   
3.    currentatom=atomsfetched(iz);
   
4. current_atom_xyz= xyz_num(iz,:);% 读入当前需要写的坐标
   
5. crt_atm_xyz_str = sprintf('%f,',current_atom_xyz);%撰写为逗号分开的string格式
   
6. crt_atm_xyz_str=crt_atm_xyz_str(1:end-1);%散掉末尾多余的逗号
   
7.   currentatom.setAttribute('XYZ', crt_atm_xyz_str);%写入对象中，atom对象直接饮用的是xsdfile这个 对象
   
8. 
   
9. end

复制代码

sobereva

我一般用gview建模，保存成gjf文件后把坐标部分拷到.mop里，算完了之后用mopac2xyz转换成xyz文件再用VMD观看。

fhh2626

sobereva 发表于 2014-12-31 05:28
我一般用gview建模，保存成gjf文件后把坐标部分拷到.mop里，算完了之后用mopac2xyz转换成xyz文件再用VMD观 ...

直接用PDBOUT关键词不就行了么。。。

sobereva

fhh2626 发表于 2014-12-31 20:54
直接用PDBOUT关键词不就行了么。。。

pdbout只能用于蛋白质（或者自己添上必要的信息），否则没法输出，不方便。
而且mopac2xyz转出来的xyz能便利地看优化轨迹，大体系的优化是mopac时下真正的用处，这需要关注优化过程。

北纬18°

fhh2626 发表于 2014-12-31 20:54
直接用PDBOUT关键词不就行了么。。。

还有这个关键，太好了！！！

北纬18°

sobereva 发表于 2015-1-1 06:54
pdbout只能用于蛋白质（或者自己添上必要的信息），否则没法输出，不方便。
而且mopac2xyz转出来的xyz ...

ubuntu 12.04.5, mopac2xyz　用不了　！

sobereva

北纬18° 发表于 2015-2-6 21:52
ubuntu 12.04.5, mopac2xyz　用不了　！

压缩包里有源代码，自行编译即可

一声叹息010

这好多个代码怎么互相调用的，能举个例子吗