diff --git a/README.md b/README.md index 9bc7afe..6a9539f 100644 --- a/README.md +++ b/README.md @@ -7,13 +7,14 @@ (2)loglammps toolkit is mainly used to extract and process data with statistic average method. (3)species toolkit can process the simulation generated species files. It sorts out the formed products with time sequence and can further refine them by molecular weight, elements and so on. Moreover, the evolution of number-average molecular weight and weight-average molecular weight can be calculated. (4)bonds_analysis toolkit can read and handle the bond order (BO) information, further mainly cope with the BO information with the following purpose: a)obtain the BO information of the specific molecule or fragment; b)give the chemical composition of all the molecules and fragments (molecular formula); c)rearrange molecular formula and trace the BO information of the interested species. -(5)lammpstrj2xyz_car_mdf_arc toolkit can output the *.xyz, *.car and *.mdf file, which can be be used for visualization by virture of Materials studio and VMD software. +(5)lammpstrj2xyz_arc_pdb toolkit can output the *.xyz, *.arc and *.pdb file, which can be be used for visualization by virture of Materials studio and VMD software, both in the form of static image and dynamic trajectory. (6)chemi_mechanism toolkit can export the products in specific trajectories and analysis the reaction path/channel. More desired functions like orientation analysis will be added in the furture. (7)structure_analysis toolkit is aimed to analysis the complex molecular structure, such as benzene ring, phenolic hydroxylic group and other specific chemical bonds or groups. However, its function is incomplete. Much work will be carried out in the furture. -3 Please cite the following papers accordingly: +3 Please cite the following papers and RMD_digging ToolKit accordingly: (1) Liu, Q.; Liu, S.; Lv, Y.; Hu, P.; Huang, Y.; Kong, M.; Li, G. Atomic-scale insight into the pyrolysis of polycarbonate by ReaxFF-based reactive molecular dynamics simulation. Fuel 2021, 287, 119484, DOI: https://doi.org/10.1016/j.fuel.2020.119484. -(2) Liu, Q.; Huang, W.; Liu, B.; Wang, P.-C.; Chen, H.-B. Gamma Radiation Chemistry of Polydimethylsiloxane Foam in Radiation-Thermal Environments: Experiments and Simulations. ACS Appl. Mat. Interfaces 2021, 13 (34), 41287-41302, DOI: https://doi.org/10.1021/acsami.1c10765. +(2) Liu, Q.; Huang, W.; Liu, B.; Wang, P.-C.; Chen, H.-B. Gamma Radiation Chemistry of Polydimethylsiloxane Foam in Radiation-Thermal Environments: Experiments and Simulations. ACS Appl. Mat. Interfaces 2021, 13 (34), 41287-41302, DOI: https://doi.org/10.1021/acsami.1c10765. +(3) Liu, Q., RMD_Digging ToolKit, https://github.com/dadaoqiuzhi/RMD_Digging. 4 Disclaimer: Using the software and the scripts in this repository is at your own risk. The software is freeware and it come without any (as in nothing, nada, niente, rien) WARRANTY. Therefore, Neither me nor the Institute of Nuclear Physics and Chemistry is responsible for the loss of data, time, bad results or anything else deriving by the use of my software or this repository. We do not make any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any software, scipts, information, apparatus, product, disclosed process, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial software, product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the People's Republic of China. The views and opinions of authors expressed herein do not necessarily state or reflect those of the People's Republic of China or the Institute of Nuclear Physics and Chemistry, and shall not be used for advertising or product endorsement purposes. \ No newline at end of file diff --git a/chemi_mechanism/.bonds_analysis_speedup.m.un~ b/chemi_mechanism/.bonds_analysis_speedup.m.un~ deleted file mode 100644 index 9ad8b3b..0000000 Binary files a/chemi_mechanism/.bonds_analysis_speedup.m.un~ and /dev/null differ diff --git a/chemi_mechanism/bonds_analysis_speedup.m~ b/chemi_mechanism/bonds_analysis_speedup.m~ deleted file mode 100644 index 304114c..0000000 --- a/chemi_mechanism/bonds_analysis_speedup.m~ +++ /dev/null @@ -1,123 +0,0 @@ -%scrit file name bonds_analysis_speedup -%purpose: -%This program is used to analyze bonds file and is more fast than -%bonds_analysis program -%version 1;2018.6.25 -disp('##################################################################################################################################') -disp('Welcome!--by Qiang Liu @Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics; Email: liubinqiang@163.com'); -disp('Repository adress of the Source code on github: https://github.com/dadaoqiuzhi/RMD_Digging'); -disp('References: 1.Fuel 287 (2021) 119484. 2.ACS Appl. Mat. Interfaces 13(34) (2021) 41287-41302. More work is coming!') -disp('##################################################################################################################################') -tartrajectory={tartrajectory(1)}; -if mod(tartrajectory{1,1},trajper)~=0 - control=0; - error('\nnonexistent trajectory, please check it!!!\n') -else - control=1; -end - -rawdata=fopen(datanamebond,'r'); -readline=0; -gap=8+atomnum; -dataline=fgetl(rawdata); -readline=readline+1; -datacell=textscan(dataline,'%s','delimiter','\n'); -datacellchar=char(datacell{1}); -datadel=strrep(datacellchar,'#',''); -datarep=strtrim(datadel); -datasplit=strsplit(datarep); -if str2num(datasplit{1,2})==tartrajectory{1} - control=0; -else - while control - i=1; - unfound=1; - while unfound - dataline=fgetl(rawdata); - readline=readline+1; - i=i+1; - if i==gap+1 - unfound=0; - break; - end - end - if mod(readline-1,gap)==0 - datacell=textscan(dataline,'%s','delimiter','\n'); - datacellchar=char(datacell{1}); - datadel=strrep(datacellchar,'#',''); - datarep=strtrim(datadel); - datasplit=strsplit(datarep); - if str2num(datasplit{1,2})==tartrajectory{1} - control=0; - end - else - disp('This is not a line with timestep, please check it!!!') - return; - end - end -end - -found=6; -while found - dataline=fgetl(rawdata); - readline=readline+1; - found=found-1; -end - -bondoutdata={};% -bondoutdata{1,1}='Timestep'; -bondoutdata{1,2}=tartrajectory{1}; -for i=3:15 - bondoutdata{1,i}=[]; -end - -line=2;atomnumcopy=atomnum; -while atomnumcopy - dataline=fgetl(rawdata); - readline=readline+1; - atomnumcopy=atomnumcopy-1; - if atomnumcopy<=0 - break; - end - datacell=textscan(dataline,'%s','delimiter','\n'); - datacellchar=char(datacell{1}); - datarep=strtrim(datacellchar); - datasplit=strsplit(datarep); - bondnumdata={}; - bondnumdata(1,1:3)=datasplit(1,1:3); - if ~strcmp(datasplit{1,3},'0') - for i=1:str2num(datasplit{1,3}) - bondnumdata(1,i+3)=datasplit(1,i+3); - end - bondnumdata(1,8)=datasplit(1,i+4); - k=i+5; - for j=1:str2num(datasplit{1,3}) - bondnumdata(1,j+8)=datasplit(1,k); - k=k+1; - end - bondnumdata(1,13:15)=datasplit(1,k:k+2); - else - bondnumdata(1,8)=datasplit(1,4); - bondnumdata(1,13)=datasplit(1,5); - bondnumdata(1,14)=datasplit(1,6); - bondnumdata(1,15)=datasplit(1,7); - end - - for kk=1:length(bondnumdata) - if isempty(bondnumdata{kk}) - bondnumdata{kk}='NaN'; - else - bondnumdata{kk}=str2num(bondnumdata{kk}); - end - end - for kk=1:length(bondnumdata) - bondoutdata{line,kk}=bondnumdata{kk}; - end - line=line+1; -end -fclose(rawdata); - -fprintf('\nbonds_analysis_speedup is successfully finished, BO information is saved inbondoutdata, search line number is recorded in readline,') - -clear atomnumcopy ans bondnumdata control datacell datacellchar datadel dataline datarep datasplit found gap i j k kk line -clear outputans unfound dataoutrow dataoutcol dataoutputrow dataoutcolchar dataoutputcol filename \ No newline at end of file diff --git a/chemi_mechanism/chemi_mechanism.m b/chemi_mechanism/chemi_mechanism.m index 4e38bb1..216766b 100644 --- a/chemi_mechanism/chemi_mechanism.m +++ b/chemi_mechanism/chemi_mechanism.m @@ -4,7 +4,7 @@ %molecular formula with the datafile of species and bonds %include modified species_analysis, species_capture, %bonds_analysis_speedup, bondorder_deepmining, lammpstrj_analysis and -%xyz_car_mdf_filemaker +%xyz_car_pdb_filemaker %version 1;2018.10.24 disp('##################################################################################################################################') disp('Welcome!--by Qiang Liu @Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics; Email: liubinqiang@163.com'); @@ -14,30 +14,26 @@ fprintf('This program is used to\n1.Export products or reactants files\n2.Analyze and export reactants-products files\n3.Analyze special groups, bonds and species (in development)\n4.Analyze where the interested species go\n') choi=input('Please select the option No.: \n'); if choi==1 || choi==2 || choi==4 - formatout=input('\nPlease select the file format No.: 1.xyz 2.arc\n'); + formatout=input('\nPlease select the file format No.: 1.xyz 2.arc 3.pdb\n'); datanamespe=input('File name of species file: \n','s'); species=input('Please input the molecular formula, element sequence should be consistent with the in.* file, multi input should be seperated by whith space: \n','s'); species=upper(species); end datanamebond=input('File name of bonds file: \n','s'); +if choi==1 || choi==2 || choi==4 + datanametrj=input('\nFile name of lammpstrj file: \n','s'); +end trajper=input('\nPlease input the output frequency of BO information and trajectory file (Positive integer, see bonds or lammpstrj file): \n'); atomnum=input('\nPlease input atom number: \n'); -fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); -if 262143>=atomnum && atomnum>32767 - fprintf('64 base number system is recommended for atom id');base=64; -elseif 32767>=atomnum && atomnum>4095 - fprintf('32 base number system is recommended for atom id');base=32; -elseif 4095>=atomnum && atomnum>999 - fprintf('16 base number system is recommended for atom id');base=16; -elseif 999>=atomnum - fprintf('10 base number system is recommended for atom id'); -elseif atomnum<0 || atomnum>32767 - error('atom number is less than 0 or larger than 32767. If larger, please check it and modify code accordingly!') -end -base=input('\nPlease select number system for atom id, positive integer and not less than the recommended value: \n'); elementsequence=input('\nPlease input atom type like C,H.O,N, seperated by white space, corresponding to 1,2,3,4...n (see *.data or in.*, \nespecially for element mapping:\n','s'); elementsequence=strtrim(elementsequence);elementsequence=strsplit(elementsequence);elementsequence=upper(elementsequence); +elemax=0; +for i=1:length(elementsequence) + if length(elementsequence{i})>elemax + elemax=length(elementsequence{i}); + end +end eleswapans=input('\nDoes there exist element mapping?y/n: \n','s'); if strcmpi(eleswapans,'y') fprintf('\nPlease input element mapping cell, eg.{''Si'',''C'';''S'',''O''},\nsingle quotes are used in practical,element in data file--expected mapping elements: \n'); @@ -46,10 +42,71 @@ else eleswap={'Nan'}; end - -if choi==1 || choi==2 || choi==4 - datanametrj=input('\nFile name of lammpstrj file: \n','s'); +if elemax==2 + if formatout==1 || formatout==2 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 262143>=atomnum && atomnum>32767 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 32767>=atomnum && atomnum>4095 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 4095>=atomnum && atomnum>999 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 999>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>32767 + error('atom number is less than 0 or larger than 262143. If larger, please check it and modify code accordingly!') + end + elseif formatout==3 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 65535>=atomnum && atomnum>16383 + fprintf('256 base number system is recommended for atom id');base=256; + elseif 16383>=atomnum && atomnum>4095 + fprintf('128 base number system is recommended for atom id');base=128; + elseif 4095>=atomnum && atomnum>1023 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 1023>=atomnum && atomnum>255 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 255>=atomnum && atomnum>99 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 99>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>65535 + error('atom number is less than 0 or larger than 65535. If larger, please check it and modify code accordingly!') + end + end +elseif elemax==1 + if formatout==1 || formatout==2 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 16777215>=atomnum && atomnum>1048575 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 1048575>=atomnum && atomnum>65535 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 65535>=atomnum && atomnum>9999 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 9999>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>16777215 + error('atom number is less than 0 or larger than 16777215. If larger, please check it and modify code accordingly!') + end + elseif formatout==3 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 262143>=atomnum && atomnum>32767 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 32767>=atomnum && atomnum>4095 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 4095>=atomnum && atomnum>999 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 999>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>262143 + error('tom number is less than 0 or larger than 262143. If larger, please check it and modify code accordingly!') + end + end +else + error('Length of atom type is NOT 1 or 2, please check it!'); end +base=input('\nPlease select number system for atom id, positive integer and not less than the recommended value: \n'); + if choi==1 maxchoi=input('\nPleasw input threshold trajectory frame value, when exceeding it different methods are suggested to limit the exportation of files: \n'); end @@ -72,12 +129,23 @@ disp('Illegal periodic boundary condition in PBCchoi, please check it!!!'); return; end - fprintf('\nScaled coordination is not recommended (use "dump_modify scale no" to avoid)') - BOXsize=input('\nDoes the coordination is scaled in the *.lammpstrj file, y/n: \n','s'); - if ~ismember(BOXsize,{'y','n'}) - error('Illegal BOXsize parameters, please check it!!!'); + elseif formatout==3 + PBCchoi=input('\nPlease input periodic boundary condition, ON/OFF: \n','s');PBCchoi=upper(PBCchoi); + if strcmp(PBCchoi,'ON') + PBCalpha=input('Periodic boundary condition, alpha, four decimal digits: \n'); + PBCbeta=input('Periodic boundary condition, beta, four decimal digits:\n'); + PBCgamma=input('Periodic boundary condition, gamma, four decimal digits:\n'); + spacegroupname=input('Point group name, eg."(P1)": \n','s');spacegroupname=upper(spacegroupname); + spacegroupname=input('\nPoint group name, eg."(P1)" for *.arc and "P 1" for *.pdb: \n','s');spacegroupname=upper(spacegroupname); + elseif strcmp(PBCchoi,'OFF') + warndlg('Nonperiodic constraint: PBC condition is not used and coordination is directly abstracted!') end end + fprintf('\nScaled coordination is not recommended (use "dump_modify scale no" to avoid)') + BOXsize=input('\nDoes the coordination is scaled in the *.lammpstrj file, y/n: \n','s');BOXsize=lower(BOXsize); + if ~ismember(BOXsize,{'y','n'}) + error('Illegal BOXsize parameters, please check it!!!'); + end if exist('outputdata','var') fprintf('\nWarning: species_analysis is not running due to the existence of outputdata in work space. \nPlease make sure this is expected species file'); @@ -304,8 +372,10 @@ datanamecar=strcat(species{1},'-',num2str(frame{1,irow}),'-',num2str(frame{2,irow}),'.xyz'); elseif formatout==2 datanamecar=strcat(species{1},'-',num2str(frame{1,irow}),'-',num2str(frame{2,irow}),'.car'); + elseif formatout==3 + datanamecar=strcat(species{1},'-',num2str(frame{1,irow}),'-',num2str(frame{2,irow}),'.pdb'); end - xyz_car_mdf_filemaker + xyz_car_pdb_filemaker seekBOinform else fclose(rawdatatrj); @@ -386,8 +456,10 @@ datanamecar=strcat(species{1},'-',num2str(frame{1,irow}),'-',num2str(frame{2,irow}),'.xyz'); elseif formatout==2 datanamecar=strcat(species{1},'-',num2str(frame{1,irow}),'-',num2str(frame{2,irow}),'.car'); + elseif formatout==3 + datanamecar=strcat(species{1},'-',num2str(frame{1,irow}),'-',num2str(frame{2,irow}),'.pdb'); end - xyz_car_mdf_filemaker + xyz_car_pdb_filemaker seekBOinform else fclose(rawdatatrj); @@ -634,32 +706,40 @@ datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','reactant','.xyz'); elseif formatout==2 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','reactant','.car'); + elseif formatout==3 + datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','reactant','.pdb'); end elseif choi==4 if formatout==1 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','product','.xyz'); elseif formatout==2 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','product','.car'); + elseif formatout==3 + datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','product','.pdb'); end end tarBOinform=tarBOinformcopy3; - xyz_car_mdf_filemaker + xyz_car_pdb_filemaker if choi==2 if formatout==1 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','product','.xyz'); elseif formatout==2 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','product','.car'); + elseif formatout==3 + datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','product','.pdb'); end elseif choi==4 if formatout==1 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','reactant','.xyz'); elseif formatout==2 datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','reactant','.car'); + elseif formatout==3 + datanamecar=strcat(species{1},'-',num2str(tartrajectory{1,1}),'-',num2str(tartrajectorycopy(1)),'-',num2str(pronum),'-','reactant','.pdb'); end end tarBOinform=tarBOinformcopy2;trjdata=trjdatacopy; - xyz_car_mdf_filemaker + xyz_car_pdb_filemaker seekBOinform diff --git a/chemi_mechanism/lammpstrj_analysis.m b/chemi_mechanism/lammpstrj_analysis.m index 3067a84..955cac9 100644 --- a/chemi_mechanism/lammpstrj_analysis.m +++ b/chemi_mechanism/lammpstrj_analysis.m @@ -89,14 +89,19 @@ for i=1:length(datasplit) boxsize(3,i)=str2num(datasplit{1,i}); end -if formatout==2 +if formatout==2 || formatout==3 if strcmp(PBCchoi,'ON') PBC='PBC=ON'; PBCa=boxsize(1,2)-boxsize(1,1); PBCb=boxsize(2,2)-boxsize(2,1); PBCc=boxsize(3,2)-boxsize(3,1); elseif strcmp(PBCchoi,'OFF') - PBC='PBC=OFF'; + if formatout==2 + PBC='PBC=OFF'; + elseif formatout==3 + PBCa=0.00;PBCb=PBCa;PBCc=PBCa; + PBCalpha=0.00;PBCbeta=PBCalpha;PBCgamma=PBCalpha; + end else disp('Illegal periodic boundary condition in PBCchoi, please check it!!!'); return; diff --git a/chemi_mechanism/xyz_car_mdf_filemaker.m b/chemi_mechanism/xyz_car_pdb_filemaker.m similarity index 53% rename from chemi_mechanism/xyz_car_mdf_filemaker.m rename to chemi_mechanism/xyz_car_pdb_filemaker.m index 364867c..a11b9da 100644 --- a/chemi_mechanism/xyz_car_mdf_filemaker.m +++ b/chemi_mechanism/xyz_car_pdb_filemaker.m @@ -11,12 +11,14 @@ disp('Repository adress of the Source code on github: https://github.com/dadaoqiuzhi/RMD_Digging'); disp('References: 1.Fuel 287 (2021) 119484. 2.ACS Appl. Mat. Interfaces 13(34) (2021) 41287-41302. More work is coming!') disp('##################################################################################################################################') -fprintf('\nxyz_car_mdf_filemaker is running, please wait...') +fprintf('\nxyz_car_pdb_filemaker is running, please wait...') if formatout==1 - title='xyz_car_mdf_filemaker Program Generated XYZ File'; + title='!xyz_car_pdb_filemaker Program Generated XYZ File'; elseif formatout==2 fileheader='!BIOSYM archive 3'; - title='xyz_car_mdf_filemaker Program Generated CAR File'; + title='!xyz_car_pdb_filemaker Program Generated CAR File'; +elseif formatout==3 + title='!xyz_car_pdb_filemaker Program Generated PDB File'; end date=datestr(now,31);date=strcat('!DATE',date); @@ -37,6 +39,16 @@ else fprintf(fid,'%s\n%s\n%-64s\n%s',fileheader,PBC,title,date); end +elseif formatout==3 + [~,tarmatchcol]=size(tarelenummatch);numofmolecule=tarmatchcol/2; + if strcmp(PBCchoi,'ON') + fprintf(fid,'%-64s\n%s\n%-9s%-9.3f%-9.3f%-8.3f%-7.2f%-7.2f%-6.2f%-14s%d',title,date,'CRYST1',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,spacegroupname,numofmolecule); + elseif strcmp(PBCchoi,'OFF') + PBCa=0.00;PBCb=PBCa;PBCc=PBCa; + PBCalpha=0.00;PBCbeta=PBCalpha;PBCgamma=PBCalpha; + [~,tarmatchcol]=size(tarelenummatch);numofmolecule=tarmatchcol/2; + fprintf(fid,'%-64s\n%s\n%-9s%-9.3f%-9.3f%-8.3f%-7.2f%-7.2f%-6.2f%-14s%d',title,date,'CRYST1',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,'None',numofmolecule); + end end element=elementsequence; numseq={}; @@ -45,13 +57,17 @@ end [tarBOrow,~]=size(tarBOinform); -readline=1; +readline=1;MOLE=1; while readline<=tarBOrow if strcmp(tarBOinform{readline,1},'#') if formatout==2 fprintf(fid,'\n%s','end'); end + if formatout==3%每个分子式有TER间隔 + fprintf(fid,'\n%3s %5d %3s %s%4d%s','TER',MOLE,'MOL','A',MOLE,'A'); + end readline=readline+1; + MOLE=MOLE+1; else elementname=charnum_match(element,numseq,tarBOinform{readline,2}); if ismember(elementname,eleswap(:,1)) @@ -67,19 +83,16 @@ tartrjdata(1,:)=trjdata(i,:); end end - if formatout==1 + if strcmpi(BOXsize,'n') xcoord=tartrjdata(3);ycoord=tartrjdata(4);zcoord=tartrjdata(5); - elseif formatout==2 - if strcmpi(BOXsize,'n') - xcoord=tartrjdata(3);ycoord=tartrjdata(4);zcoord=tartrjdata(5); - elseif strcmpi(BOXsize,'y') - xcoord=boxsize(1,1)+tartrjdata(3)*PBCa; - ycoord=boxsize(2,1)+tartrjdata(4)*PBCb; - zcoord=boxsize(3,1)+tartrjdata(5)*PBCc; - else - disp('Illegal BOXsize parameters!!!') - end - + elseif strcmpi(BOXsize,'y') + xcoord=boxsize(1,1)+tartrjdata(3)*PBCa; + ycoord=boxsize(2,1)+tartrjdata(4)*PBCb; + zcoord=boxsize(3,1)+tartrjdata(5)*PBCc; + else + disp('Illegal BOXsize parameters!!!') + end + if formatout==2 ff=forcefield(element,numseq,tarBOinform{readline,2}); if ismember(upper(ff),eleswap(:,1)) [~,lib]=ismember(upper(ff),eleswap(:,1)); @@ -91,16 +104,43 @@ fprintf(fid,'\n%-5s %14.09f %14.09f %14.09f',atomname,xcoord,ycoord,zcoord); elseif formatout==2 fprintf(fid,'\n%-5s %14.09f %14.09f %14.09f %-12s%-7s %-2s %6.3f',atomname,xcoord,ycoord,zcoord,residueseqname,ff,elementname,charge); + elseif formatout==3 + atomNO=tartrjdata(1); + charge=tarBOinform{readline,15}; + fprintf(fid,'\n%4s %5d %4s %3s %4d %8.03f%8.03f%8.03f%6.02f%6.02f %2s%2.01f','ATOM',atomNO,atomname,'MOL',MOLE,xcoord,ycoord,zcoord,1.00,0.00,elementname,charge); end readline=readline+1; - atomcount=atomcount+1; + if formatout==1 + atomcount=atomcount+1; + end end end if formatout==2 - fprintf(fid,'\n%s\n','end'); + fprintf(fid,'\n%s','end'); +end +if formatout==3 + %tarBOinform中查找连接信息 + trjreadline=1; + while trjreadline<=tarBOrow + if tarBOinform{trjreadline,3}==0 + fprintf(fid,'\n%6s%5d','CONECT',tarBOinform{trjreadline,1}); + elseif tarBOinform{trjreadline,3}==1 + fprintf(fid,'\n%6s%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4}); + elseif tarBOinform{trjreadline,3}==2 + fprintf(fid,'\n%6s%5d%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4},tarBOinform{trjreadline,5}); + elseif tarBOinform{trjreadline,3}==3 + fprintf(fid,'\n%6s%5d%5d%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4},tarBOinform{trjreadline,5},tarBOinform{trjreadline,6}); + elseif tarBOinform{trjreadline,3}==4 + fprintf(fid,'\n%6s%5d%5d%5d%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4},tarBOinform{trjreadline,5},tarBOinform{trjreadline,6},tarBOinform{trjreadline,7}); + elseif ~isnumeric(tarBOinform{trjreadline,3}) && ~strcmpi(tarBOinform{trjreadline,3},'#') + error('Number of chemical bond is NOT between 0 and 4, please check it!!!'); + end + trjreadline=trjreadline+1; + end + fprintf(fid,'\n%s','END'); end fclose(fid); -fprintf('\nxyz_car_mdf_filemaker is successfully finished.') +fprintf('\nxyz_car_pdb_filemaker is successfully finished.') clear ans atomname charge element elementname ff fid numseq clear tarrow xcoord ycoord zcoord topo topocol atomcount diff --git a/lammpstrj2xyz_car_mdf_arc/SysConvert.m b/lammpstrj2xyz_arc_pdb/SysConvert.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/SysConvert.m rename to lammpstrj2xyz_arc_pdb/SysConvert.m diff --git a/lammpstrj2xyz_car_mdf_arc/atomnummolecule_strcat.m b/lammpstrj2xyz_arc_pdb/atomnummolecule_strcat.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/atomnummolecule_strcat.m rename to lammpstrj2xyz_arc_pdb/atomnummolecule_strcat.m diff --git a/lammpstrj2xyz_car_mdf_arc/bondorder_deepmining.m b/lammpstrj2xyz_arc_pdb/bondorder_deepmining.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/bondorder_deepmining.m rename to lammpstrj2xyz_arc_pdb/bondorder_deepmining.m diff --git a/lammpstrj2xyz_car_mdf_arc/bonds_analysis_speedup.m b/lammpstrj2xyz_arc_pdb/bonds_analysis_speedup.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/bonds_analysis_speedup.m rename to lammpstrj2xyz_arc_pdb/bonds_analysis_speedup.m diff --git a/lammpstrj2xyz_car_mdf_arc/cellrowcol_del.m b/lammpstrj2xyz_arc_pdb/cellrowcol_del.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/cellrowcol_del.m rename to lammpstrj2xyz_arc_pdb/cellrowcol_del.m diff --git a/lammpstrj2xyz_car_mdf_arc/charnum_match.m b/lammpstrj2xyz_arc_pdb/charnum_match.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/charnum_match.m rename to lammpstrj2xyz_arc_pdb/charnum_match.m diff --git a/lammpstrj2xyz_car_mdf_arc/eleme_molecule.m b/lammpstrj2xyz_arc_pdb/eleme_molecule.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/eleme_molecule.m rename to lammpstrj2xyz_arc_pdb/eleme_molecule.m diff --git a/lammpstrj2xyz_car_mdf_arc/forcefield.m b/lammpstrj2xyz_arc_pdb/forcefield.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/forcefield.m rename to lammpstrj2xyz_arc_pdb/forcefield.m diff --git a/lammpstrj2xyz_car_mdf_arc/lammpstrj_analysis.m b/lammpstrj2xyz_arc_pdb/lammpstrj_analysis.m similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/lammpstrj_analysis.m rename to lammpstrj2xyz_arc_pdb/lammpstrj_analysis.m diff --git a/lammpstrj2xyz_car_mdf_arc/output_mydata.xlsx b/lammpstrj2xyz_arc_pdb/output_mydata.xlsx similarity index 100% rename from lammpstrj2xyz_car_mdf_arc/output_mydata.xlsx rename to lammpstrj2xyz_arc_pdb/output_mydata.xlsx diff --git a/lammpstrj2xyz_car_mdf_arc/xyz_arc_filemaker_speedupMOLE.m b/lammpstrj2xyz_arc_pdb/xyz_arc_filemaker_speedupMOLE.m similarity index 97% rename from lammpstrj2xyz_car_mdf_arc/xyz_arc_filemaker_speedupMOLE.m rename to lammpstrj2xyz_arc_pdb/xyz_arc_filemaker_speedupMOLE.m index cf83ff8..4d62fde 100644 --- a/lammpstrj2xyz_car_mdf_arc/xyz_arc_filemaker_speedupMOLE.m +++ b/lammpstrj2xyz_arc_pdb/xyz_arc_filemaker_speedupMOLE.m @@ -114,9 +114,9 @@ error('Illegal BOXsize parameters, please check it!!!'); end if formatout==1 - title='xyz_arc_filemaker_speedupMOLE Program Generated XYZ File'; + title='!xyz_arc_filemaker_speedupMOLE Program Generated XYZ File'; elseif formatout==2 - title='xyz_arc_filemaker_speedupMOLE Program Generated ARC File'; + title='!xyz_arc_filemaker_speedupMOLE Program Generated ARC File'; end date=datestr(now,31);date=strcat('!DATE',date); @@ -132,7 +132,7 @@ disp('Periodic boundary condition, eg."PBC 33.4531 33.4531 33.4531 90.0000 90.0000 90.0000 (P1)"'); PBCgamma=input('Periodic boundary condition, gamma, four decimal digits:\n'); disp('Periodic boundary condition, eg."PBC 33.4531 33.4531 33.4531 90.0000 90.0000 90.0000 (P1)"'); - spacegroupname=input('\nPoint group name, eg."(P1)": \n,'s');spacegroupname=upper(spacegroupname); + spacegroupname=input('\nPoint group name, eg."(P1)": \n','s');spacegroupname=upper(spacegroupname); fprintf(fid,'%s\n%s',fileheader,PBC); elseif strcmp(PBCchoi,'OFF') PBC='PBC=OFF'; @@ -314,7 +314,6 @@ continue end fprintf('\nStep2:Group %d trajectory %d is successfully processed by bondorder_deepmining, and tarBOinform is generated,\ncontinue running lammpstrj_analysis, please wait...\n',ii,trjcollection{1,ii}); -fprintf('\nStep1:Group %d trajectory %d is successfully processed by bonds_analysis_speedup and bondnumdata is generated,\n continue running bondorder_deepmining program, please wait...\n',ii,trjcollection{1,ii}); @@ -478,7 +477,7 @@ end end if formatout==2 - fprintf(fid,'\n%s\n','end'); + fprintf(fid,'\n%s','end'); end diff --git a/lammpstrj2xyz_arc_pdb/xyz_arc_pdb_filemaker.m b/lammpstrj2xyz_arc_pdb/xyz_arc_pdb_filemaker.m new file mode 100644 index 0000000..d3bb178 --- /dev/null +++ b/lammpstrj2xyz_arc_pdb/xyz_arc_pdb_filemaker.m @@ -0,0 +1,629 @@ +%scrit file name xyz_arc_pdb_filemaker +%purpose: +%This program is used to generate .arc file which containing a sequence of structures for bonds and *.lammpstrj file of REAXC +%include bondoutdata generated by bonds_analysis program in SIB +%mode or produced by bonds_analysis_speedup program (recommend) +%include tarBOinform and tarelenummatch file generated by bondorder_deepmining program +%inclunde trjdata file generated by lammpstrj_analysis program +%version 1;2021.10.2 +disp('##################################################################################################################################') +disp('Welcome!--by Qiang Liu @Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics; Email: liubinqiang@163.com'); +disp('Repository adress of the Source code on github: https://github.com/dadaoqiuzhi/RMD_Digging'); +disp('References: 1.Fuel 287 (2021) 119484. 2.ACS Appl. Mat. Interfaces 13(34) (2021) 41287-41302. More work is coming!') +disp('##################################################################################################################################') +fprintf('\nThis program will turn BO information in tarBOinform and trjdata into.xyz,.arc or .pdb file, which can be visualized by VMD or Materials studio.\nCaution: Write in attach mode!!!\n') +formatout=input('\nPlease select the option No. of file format: 1.xyz 2.arc 3.pdb\n'); +if formatout==1 + dataname=input('\nPlease give the file name *.xyz, e.g."PCdeg.xyz"\n','s'); +elseif formatout==2 + dataname=input('\nPlease give the file name *.arc, e.g."PCdeg.arc"\n','s'); +elseif formatout==3 + dataname=input('\nPlease give the file name *.pdb, e.g."PCdeg.pdb"\n','s'); +else + error('Illegal file format No., please check it!!!') +end +bonddataname=input('\nFilename of bond file: \n','s'); +bondper=input('\nPlease input the output frequency of BO information (Positive integer, can be confirmed in bond file): \n'); +trjdataname=input('\nFilename of *.lammpstrj file: \n','s'); +trjper=input('\nPlease input the output frequency of BO information and trajectory file (Positive integer, see bonds or lammpstrj file):\n'); +atomnum=input('Please input atom number: \n'); + +elementsequence=input('\nPlease input atom type like C,H.O,N, seperated by white space, corresponding to 1,2,3,4...n (see *.data or in.*, \nespecially for element mapping:\n','s'); +element=upper(elementsequence);element=strtrim(element);element=strsplit(element); +numseq={};elemax=0; +for i=1:length(element) + numseq{1,i}=i; + if length(element{i})>elemax + elemax=length(element{i}); + end +end +eleswapans=input('\nDoes there exist element mapping?y/n: \n','s'); +if strcmpi(eleswapans,'y') + fprintf('\nPlease input element mapping cell, eg.{''Si'',''C'';''S'',''O''},\nsingle quotes are used in practical,element in data file--expected mapping elements: \n'); + eleswap=input(''); + eleswap=upper(eleswap); +else + eleswap={'Nan'}; +end + +if elemax==2 + if formatout==1 || formatout==2 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 262143>=atomnum && atomnum>32767 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 32767>=atomnum && atomnum>4095 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 4095>=atomnum && atomnum>999 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 999>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>32767 + error('atom number is less than 0 or larger than 262143. If larger, please check it and modify code accordingly!') + end + elseif formatout==3 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 65535>=atomnum && atomnum>16383 + fprintf('256 base number system is recommended for atom id');base=256; + elseif 16383>=atomnum && atomnum>4095 + fprintf('128 base number system is recommended for atom id');base=128; + elseif 4095>=atomnum && atomnum>1023 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 1023>=atomnum && atomnum>255 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 255>=atomnum && atomnum>99 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 99>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>65535 + error('atom number is less than 0 or larger than 65535. If larger, please check it and modify code accordingly!') + end + end +elseif elemax==1 + if formatout==1 || formatout==2 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 16777215>=atomnum && atomnum>1048575 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 1048575>=atomnum && atomnum>65535 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 65535>=atomnum && atomnum>9999 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 9999>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>16777215 + error('atom number is less than 0 or larger than 16777215. If larger, please check it and modify code accordingly!') + end + elseif formatout==3 + fprintf('\nDifferent number system is adopted according to the atom number (ASCII)\n'); + if 262143>=atomnum && atomnum>32767 + fprintf('64 base number system is recommended for atom id');base=64; + elseif 32767>=atomnum && atomnum>4095 + fprintf('32 base number system is recommended for atom id');base=32; + elseif 4095>=atomnum && atomnum>999 + fprintf('16 base number system is recommended for atom id');base=16; + elseif 999>=atomnum + fprintf('10 base number system is recommended for atom id'); + elseif atomnum<0 || atomnum>262143 + error('tom number is less than 0 or larger than 262143. If larger, please check it and modify code accordingly!') + end + end +else + error('Length of atom type is NOT 1 or 2, please check it!'); +end +base=input('\nPlease select number system for atom id, positive integer and not less than the recommended value: \n'); +trjatomnum=atomnum; +trjcollection=input('\nPlease select the export method No.: \n1.Manually specify multi-trajectories \n2.Monotonically increasing frame(s) in arithmetic sequence(closed interval) \n3.Both method 1 and2 (1 first, 2 last)\n'); +if trjcollection==1 + trjcollection=input('\nTrajectories in ascending order, seperated by white space\n','s'); + trjcollection=strtrim(trjcollection);trjcollection=strsplit(trjcollection); + for i=1:length(trjcollection) + trjcollection{1,i}=str2num(trjcollection{1,i}); + end + +elseif trjcollection==2 + trjcollection=input('\nMinimum and maximum trajectories, and equal difference: \n','s'); + trjcollection=strtrim(trjcollection);trjcollection=strsplit(trjcollection); + for i=1:length(trjcollection) + trjcollection{1,i}=str2num(trjcollection{1,i}); + end + trjmin=trjcollection{1,1};trjmax=trjcollection{1,2};trjstep=trjcollection{1,3}; + if trjmax>trjmin + trjmod=mod(trjmax-trjmin,trjstep);trjnnum=(trjmax-trjmin-trjmod)/trjstep; + else + error('Illegal trajectory, please check it!!!'); + end + for i=1:trjnnum + trjcollection{1,1+i}=trjmin+trjstep*i; + end + if trjmod~=0 + trjcollection{1,2+trjnnum}=trjmax; + end + +elseif trjcollection==3 + trjcollection=input('\nManually specification first and then Minimum and maximum trajectories, and equal difference, seperated by white space: \n','s'); + trjcollection=strtrim(trjcollection);trjcollection=strsplit(trjcollection); + lengthtrj=length(trjcollection);trjcollect={}; + for i=1:lengthtrj-3 + trjcollect{1,i}=str2num(trjcollection{1,i}); + end + trjmin=str2num(trjcollection{1,lengthtrj-2}); + trjmax=str2num(trjcollection{1,lengthtrj-1}); + trjstep=str2num(trjcollection{1,lengthtrj}); + trjcollect{1,lengthtrj-2}=str2num(trjcollection{1,lengthtrj-2}); + if trjmax>trjmin + trjmod=mod(trjmax-trjmin,trjstep);trjnnum=(trjmax-trjmin-trjmod)/trjstep; + else + error('Illegal trajectory, please check it!!!'); + end + for k=1:trjnnum + trjcollect{1,lengthtrj-2+k}=trjmin+trjstep*k; + end + if trjmod~=0 + trjcollect{1,lengthtrj-1+trjnnum}=trjmax; + end + trjcollection=trjcollect; +else + error('Illegal method No., please check it!!!'); +end + +trjcollection=cell2mat(trjcollection);trjcollection=sort(trjcollection,2);trjlength=length(trjcollection); +trjones=ones(1,trjlength);trjcollection=mat2cell(trjcollection,1,trjones); + +for i=1:length(trjcollection) + if mod(trjcollection{1,i},bondper)~=0 + error('Nonexistent trajectory in bond file, please check it!!!') + end + if mod(trjcollection{1,i},trjper)~=0 + error('Nonexistent trajectory in *.lammpstrj, please check it!!!') + end +end + +if formatout==2 + fileheader='!BIOSYM archive 3'; + PBCchoi=input('\nPlease input periodic boundary condition, ON/OFF: \n','s');PBCchoi=upper(PBCchoi); +end +if formatout==3 + PBCchoi=input('\nPlease input periodic boundary condition, ON/OFF: \n','s');PBCchoi=upper(PBCchoi); + if strcmpi(PBCchoi,'OFF') + warndlg('Nonperiodic constraint: PBC condition is not used and coordination is directly abstracted!') + end +end +BOXsize=input('\nDoes the coordination is scaled in the *.lammpstrj file, y/n: \n','s');BOXsize=lower(BOXsize); +if ~ismember(BOXsize,{'y','n'}) + error('Illegal BOXsize parameters, please check it!!!'); +end +if formatout==1 + title='!xyz_arc_pdb_filemaker Program Generated XYZ File'; +elseif formatout==2 + title='!xyz_arc_pdb_filemaker Program Generated ARC File'; +elseif formatout==3 + title='!xyz_arc_pdb_filemaker Program Generated PDB File'; +end +date=datestr(now,31);date=strcat('!DATE',date); + +fid=fopen(dataname,'at'); + +if formatout==2 || formatout==3 + if strcmp(PBCchoi,'ON') + PBC='PBC=ON'; + PBCalpha=input('Periodic boundary condition, alpha, four decimal digits: \n'); + PBCbeta=input('Periodic boundary condition, beta, four decimal digits:\n'); + PBCgamma=input('Periodic boundary condition, gamma, four decimal digits:\n'); + spacegroupname=input('\nPoint group name, eg."(P1)" for *.arc and "P 1" for *.pdb: \n','s');spacegroupname=upper(spacegroupname); + if formatout==2 + fprintf(fid,'%s\n%s',fileheader,PBC); + end + elseif strcmp(PBCchoi,'OFF') + PBC='PBC=OFF'; + if formatout==2 + fprintf(fid,'%s\n%s',fileheader,PBC); + end + else + disp('Illegal periodic boundary condition in PBCchoi, please check it!!!'); + return; + end +end + +fprintf('\nxyz_arc_pdb_filemaker is running, p;ease wait...\n') + +ii=1; +rawdata=fopen(bonddataname,'r'); +dataline=fgetl(rawdata);datacell=textscan(dataline,'%s','delimiter','\n'); +datacellchar=char(datacell{1});datadel=strrep(datacellchar,'#',''); +datarep=strtrim(datadel); + +while ii<=length(trjcollection) + fprintf('\nxyz_arc_pdb_filemaker is searching for Group %d trajectory: %d,%d trajectories in total\n',ii,trjcollection{1,ii},trjlength); + control=1;atomnum=trjatomnum; + gap=8+atomnum; + readline=1; + if ii>=2 + dataline=fgetl(rawdata); + dataline=fgetl(rawdata); + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1});datadel=strrep(datacellchar,'#',''); + datarep=strtrim(datadel); + end + datasplit=strsplit(datarep); + if str2num(datasplit{1,2})==trjcollection{1,ii} + control=0; + else + while control + dataline=textscan(rawdata,'%q',1,'headerlines',gap-1,'delimiter','\n');%ָ + readline=readline+gap; + if mod(readline-1,gap)==0 + datacell=dataline; + datacellchar=char(datacell{1}); + datadel=strrep(datacellchar,'#',''); + datarep=strtrim(datadel); + datasplit=strsplit(datarep); + if str2num(datasplit{1,2})==trjcollection{1,ii} + control=0; + end + else + error('bonds_analysis_speedup found a line does not belong to a timestep line, please check it!!!') + end + end + end + dataline=textscan(rawdata,'%q',1,'headerlines',5,'delimiter','\n'); + bondoutdata={}; + bondoutdata{1,1}='Timestep'; + bondoutdata{1,2}=trjcollection{1,ii}; + for i=3:15 + bondoutdata{1,i}=[]; + end + line=2; + while atomnum + dataline=fgetl(rawdata); + readline=readline+1; + atomnum=atomnum-1; + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + datasplit=strsplit(datarep); + bondnumdata={}; + bondnumdata(1,1:3)=datasplit(1,1:3); + if str2num(datasplit{1,3})~=0 + for i=1:str2num(datasplit{1,3}) + bondnumdata(1,i+3)=datasplit(1,i+3); + end + else + i=0; + end + bondnumdata(1,8)=datasplit(1,i+4); + k=i+5; + for j=1:str2num(datasplit{1,3}) + bondnumdata(1,j+8)=datasplit(1,k); + k=k+1; + end + bondnumdata(1,13:15)=datasplit(1,k:k+2); + for kk=1:length(bondnumdata) + if isempty(bondnumdata{kk}) + bondnumdata{kk}='NaN'; + else + bondnumdata{kk}=str2num(bondnumdata{kk}); + end + end + for kk=1:length(bondnumdata) + bondoutdata{line,kk}=bondnumdata{kk}; + end + line=line+1; + end + fprintf('\nStep1:Group %d trajectory %d is successfully processed by bonds_analysis_speedup and bondnumdata is generated,\n continue running bondorder_deepmining program, please wait...\n',ii,trjcollection{1,ii}); + + + separator={'#','#','#','#','#','#','#','#','#','#','#','#','#','#','#'}; + [row,~]=size(bondoutdata);bondoutdata(row+1,:)=separator(1,:);tarbondnum=row-1; + tarelenummatch={};tarBOinform={};lineofelenum=1;tartrjnum=1; + while ~ischar(bondoutdata{tartrjnum+1,1}) + datapython={};datapython{1,1}=bondoutdata{tartrjnum+1,1}; + for i=1:bondoutdata{tartrjnum+1,3} + datapython{1,i+1}=bondoutdata{tartrjnum+1,3+i}; + end + elenummatch={};elenummatch(:,1)=element'; + for i=1:length(element) + elenummatch{i,2}=0; + end + [elenumrow,~]=size(elenummatch); + elementname=charnum_match(element,numseq,bondoutdata{tartrjnum+1,2}); + if ismember(elementname,eleswap(:,1)) + [~,lib]=ismember(elementname,eleswap(:,1)); + elementname=eleswap{lib,2}; + end + elenummatch=eleme_molecule(elenummatch,elementname); + BOinform={};BOinform(1,:)=bondoutdata(tartrjnum+1,:); + lineofbo=2; + bondoutdata{tartrjnum+1,1}='NaN'; + datapython{1,1}='NaN'; + bondoutdata=cellrowcol_del(bondoutdata,'delrow','NaN'); + tarbondnum=tarbondnum-1; + datapython=cellrowcol_del(datapython,'delcol','NaN'); + while ~isempty(datapython) && ~ischar(bondoutdata{tartrjnum+1,1}) + alter=datapython{1,1};kk=0; + for i=tartrjnum+1:tartrjnum+tarbondnum + if alter==datapython{1,1} + kk=kk+1; + if datapython{1,1}==bondoutdata{i,1} + j=length(datapython)+1; + for k=1:bondoutdata{i,3} + datapython{1,j}=bondoutdata{i,3+k}; + j=j+1; + end + datapython{1,1}='NaN'; + datapython=cellrowcol_del(datapython,'delcol','NaN'); + elementname=charnum_match(element,numseq,bondoutdata{i,2}); + if ismember(elementname,eleswap(:,1)) + [~,lib]=ismember(elementname,eleswap(:,1)); + elementname=eleswap{lib,2}; + end + elenummatch=eleme_molecule(elenummatch,elementname); + BOinform(lineofbo,:)=bondoutdata(i,:); + lineofbo=lineofbo+1; + bondoutdata{i,1}='NaN'; + bondoutdata=cellrowcol_del(bondoutdata,'delrow','NaN'); + tarbondnum=tarbondnum-1; + break + end + if alter==datapython{1,1} && kk==tarbondnum + datapython{1,1}='NaN'; + datapython=cellrowcol_del(datapython,'delcol','NaN'); + break + end + end + end + end + tarelenummatch(1:elenumrow,2*lineofelenum-1:2*lineofelenum)=elenummatch(1:elenumrow,1:2); + lineofelenum=lineofelenum+1;[BOrow,~]=size(BOinform); + BOinform(BOrow+1,:)=separator(1,:); + [rowtarBO,~]=size(tarBOinform); + tarBOinform(rowtarBO+1:rowtarBO+BOrow+1,:)=BOinform(1:BOrow+1,:); + continue + end%bondorder_deepmining + fprintf('\nStep2:Group %d trajectory %d is successfully processed by bondorder_deepmining, and tarBOinform is generated,\ncontinue running lammpstrj_analysis, please wait...\n',ii,trjcollection{1,ii}); + + + + trjreadline=0;control=1;atomnum=trjatomnum; + gap=9+atomnum; + trjrawdata=fopen(trjdataname,'r'); + dataline=fgetl(trjrawdata); + readline=readline+1; + trjreadline=trjreadline+1; + dataline=fgetl(trjrawdata); + readline=readline+1; + trjreadline=trjreadline+1; + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + if str2num(datarep)==trjcollection{1,ii} + control=0; + else + while control + dataline=textscan(trjrawdata,'%q',1,'headerlines',gap-1,'delimiter','\n'); + trjreadline=trjreadline+gap; + if mod(trjreadline-2,gap)==0 + datacell=dataline; + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + if str2num(datarep)==trjcollection{1,ii} + control=0; + end + else + error('lammpstrj_analysis found a line does not belong to a timestep line, please check it!!!'); + end + end + end + + dataline=textscan(trjrawdata,'%q',1,'headerlines',2,'delimiter','\n'); + boxsize=[]; + dataline=fgetl(trjrawdata); + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + datasplit=strsplit(datarep); + for i=1:length(datasplit) + boxsize(1,i)=str2num(datasplit{1,i}); + end + dataline=fgetl(trjrawdata); + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + datasplit=strsplit(datarep); + for i=1:length(datasplit) + boxsize(2,i)=str2num(datasplit{1,i}); + end + dataline=fgetl(trjrawdata); + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + datasplit=strsplit(datarep); + for i=1:length(datasplit) + boxsize(3,i)=str2num(datasplit{1,i}); + end + if formatout==1 + if ii==1 + fprintf(fid,'%d\n%s %s',atomnum,title,date); + elseif ii>1 + fprintf(fid,'\n%d\n%s %s',atomnum,title,date); + end + if strcmpi(BOXsize,'y') + PBCa=boxsize(1,2)-boxsize(1,1); + PBCb=boxsize(2,2)-boxsize(2,1); + PBCc=boxsize(3,2)-boxsize(3,1); + end + elseif formatout==2 + + if strcmp(PBCchoi,'ON') + PBC='PBC=ON'; + PBCa=boxsize(1,2)-boxsize(1,1); + PBCb=boxsize(2,2)-boxsize(2,1); + PBCc=boxsize(3,2)-boxsize(3,1); + fprintf(fid,'\n%-64s\n%s\n%3s%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f %-16s',title,date,'PBC',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,spacegroupname); + elseif strcmp(PBCchoi,'OFF') + PBC='PBC=OFF'; + fprintf(fid,'\n%-64s\n%s\',title,date); + else + disp('Illegal periodic boundary condition in PBCchoi, please check it!!!'); + return; + end + elseif formatout==3 + if strcmp(PBCchoi,'ON') + PBCa=boxsize(1,2)-boxsize(1,1); + PBCb=boxsize(2,2)-boxsize(2,1); + PBCc=boxsize(3,2)-boxsize(3,1); + [~,tarmatchcol]=size(tarelenummatch);numofmolecule=tarmatchcol/2; + if ii==1 + fprintf(fid,'%-64s\n%s\n%-9s%-9.3f%-9.3f%-8.3f%-7.2f%-7.2f%-6.2f%-14s%d',title,date,'CRYST1',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,spacegroupname,numofmolecule); + else + fprintf(fid,'\n%-64s\n%s\n%-9s%-9.3f%-9.3f%-8.3f%-7.2f%-7.2f%-6.2f%-14s%d',title,date,'CRYST1',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,spacegroupname,numofmolecule); + end + elseif strcmp(PBCchoi,'OFF') + PBCa=0.00;PBCb=PBCa;PBCc=PBCa; + PBCalpha=0.00;PBCbeta=PBCalpha;PBCgamma=PBCalpha; + [~,tarmatchcol]=size(tarelenummatch);numofmolecule=tarmatchcol/2; + if ii==1 + fprintf(fid,'%-64s\n%s\n%-9s%-9.3f%-9.3f%-8.3f%-7.2f%-7.2f%-6.2f%-14s%d',title,date,'CRYST1',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,'None',numofmolecule); + else + fprintf(fid,'\n%-64s\n%s\n%-9s%-9.3f%-9.3f%-8.3f%-7.2f%-7.2f%-6.2f%-14s%d',title,date,'CRYST1',PBCa,PBCb,PBCc,PBCalpha,PBCbeta,PBCgamma,'None',numofmolecule); + end + end + end + + dataline=textscan(trjrawdata,'%q',1,'headerlines',0,'delimiter','\n'); + trjdata=[];line=1; + while atomnum%ȡϢ + dataline=fgetl(trjrawdata); + readline=readline+1; + trjreadline=trjreadline+1; + atomnum=atomnum-1; + datacell=textscan(dataline,'%s','delimiter','\n'); + datacellchar=char(datacell{1}); + datarep=strtrim(datacellchar); + datasplit=strsplit(datarep); + for i=1:length(datasplit) + trjdata(line,i)=str2num(datasplit{1,i}); + end + line=line+1; + end + if formatout==1 + fprintf('\nStep3:Group %d trajectory %d is successfully processed b lammpstrj_analysis and trjdata is generated, \ncontinue running xyz_arc_filemaker_speedupMOLE to generate *.xyz fileplease waits...\n',ii,trjcollection{1,ii}) + elseif formatout==2 + fprintf('\nStep3:Group %d trajectory %d is successfully processed b lammpstrj_analysis and trjdata is generated, \ncontinue running xyz_arc_filemaker_speedupMOLE to generate *.arc fileplease waits...\n',ii,trjcollection{1,ii}) + end + + [~,tarmatchcol]=size(tarelenummatch);numofmolecule=tarmatchcol/2; + [tarBOrow,~]=size(tarBOinform);line=1;trjreadline=1; + MOLE=1; + while trjreadline<=tarBOrow + if strcmp(tarBOinform{trjreadline,1},'#') + if formatout==2 + fprintf(fid,'\n%s','end'); + end + if formatout==3 + fprintf(fid,'\n%3s %5d %3s %s%4d%s','TER',MOLE,'MOL','A',MOLE,'A'); + end + trjreadline=trjreadline+1; + readline=readline+1; + MOLE=MOLE+1; + else + elementname=charnum_match(element,numseq,tarBOinform{trjreadline,2}); + if ismember(elementname,eleswap(:,1)) + [~,lib]=ismember(elementname,eleswap(:,1)); + elementname=eleswap{lib,2}; + end + atomid_conv=SysConvert(tarBOinform{trjreadline,1},base); + atomname=strcat(elementname,atomid_conv); + + [trjrow,~]=size(trjdata);tartrjdata=[]; + for i=1:trjrow + if trjdata(i,1)==tarBOinform{trjreadline,1} + tartrjdata(1,:)=trjdata(i,:); + end + end + if strcmpi(BOXsize,'n') + xcoord=tartrjdata(3);ycoord=tartrjdata(4);zcoord=tartrjdata(5); + elseif strcmpi(BOXsize,'y') + xcoord=boxsize(1,1)+tartrjdata(3)*PBCa; + ycoord=boxsize(2,1)+tartrjdata(4)*PBCb; + zcoord=boxsize(3,1)+tartrjdata(5)*PBCc; + else + error('Illegal BOXsize parameters, please check it!!!'); + end + + if formatout==1 + fprintf(fid,'\n%-5s %14.09f %14.09f %14.09f',atomname,xcoord,ycoord,zcoord); + elseif formatout==2 + residueseqname= strcat('MOLE',' #',num2str(MOLE)); + residueseqname=strrep(residueseqname,'#',''); + ff=forcefield(element,numseq,tarBOinform{trjreadline,2}); + if ismember(upper(ff),eleswap(:,1)) + [~,lib]=ismember(upper(ff),eleswap(:,1)); + ff=lower(eleswap{lib,2}); + end + charge=tarBOinform{trjreadline,15}; + + fprintf(fid,'\n%-5s %14.09f %14.09f %14.09f %-12s%-7s %-2s %6.3f',atomname,xcoord,ycoord,zcoord,residueseqname,ff,elementname,charge); + elseif formatout==3 + atomNO=tartrjdata(1); + charge=tarBOinform{trjreadline,15}; + fprintf(fid,'\n%4s %5d %4s %3s %4d %8.03f%8.03f%8.03f%6.02f%6.02f %2s%2.01f','ATOM',atomNO,atomname,'MOL',MOLE,xcoord,ycoord,zcoord,1.00,0.00,elementname,charge); + end + trjreadline=trjreadline+1; + readline=readline+1; + end + end + + if formatout==2 + fprintf(fid,'\n%s','end'); + end + + if formatout==3 + + trjreadline=1; + while trjreadline<=tarBOrow + if tarBOinform{trjreadline,3}==0 + fprintf(fid,'\n%6s%5d','CONECT',tarBOinform{trjreadline,1}); + elseif tarBOinform{trjreadline,3}==1 + fprintf(fid,'\n%6s%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4}); + elseif tarBOinform{trjreadline,3}==2 + fprintf(fid,'\n%6s%5d%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4},tarBOinform{trjreadline,5}); + elseif tarBOinform{trjreadline,3}==3 + fprintf(fid,'\n%6s%5d%5d%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4},tarBOinform{trjreadline,5},tarBOinform{trjreadline,6}); + elseif tarBOinform{trjreadline,3}==4 + fprintf(fid,'\n%6s%5d%5d%5d%5d%5d','CONECT',tarBOinform{trjreadline,1},tarBOinform{trjreadline,4},tarBOinform{trjreadline,5},tarBOinform{trjreadline,6},tarBOinform{trjreadline,7}); + elseif ~isnumeric(tarBOinform{trjreadline,3}) && ~strcmpi(tarBOinform{trjreadline,3},'#') + error('ѧĿ0-4֮䣬飡'); + end + trjreadline=trjreadline+1; + end + fprintf(fid,'\n%s','END'); + end + + if ii==length(trjcollection) + break; + end + ii=ii+1; +end +fclose(rawdata); +fclose(trjrawdata); +fclose(fid); +fprintf('\nxyz_arc_pdb_filemaker is successfully finished.\n\n') +if formatout==1 + fprintf('\n*.xyz file is generated: %s\n',dataname) +elseif formatout==2 + fprintf('\n*.arc file is generated: %s\n',dataname) +elseif formatout==3 + fprintf('\n*.pdb file is generated: %s\n',dataname) +end +msgbox('Donexyz_arc_pdb_filemaker is successfully finished.'); + +clear ans atomnum bondnumdata control datacell datacellchar datadel dataline dataname datarep datasplit found gap i j k kk line +clear outputans rawdata tartrajectory trajper unfound dataoutrow dataoutcol dataoutputrow dataoutcolchar dataoutputcol filename +clear alter bondrownum BOrow col datapython element elementname elementsequence elementsequence readline atomNO elemax +clear elenummatch elenumrow i j k kk lineofbo lineofelenum numseq row rowtarBO separator speciestrjnum spacegroupname +clear tarbondnum tartrjnum trajectorynum bonddataname bondper ii PBCa PBCalpha PBCb PBCbeta PBCc PBCgamma tarelenummatch +clear ans atomname BOXsize charge dataname date element elementname elementsequence ff fid fileheader numseq PBC PBCchoi PBCcond +clear residuename residueseqname tarrow title xcoord xhi xlo xlength ycoord yhi ylo ylength zcoord zhi zlo zlength date topo topocol +clear molecule comatomname mdfans symmetry groupname connect row connectivity i j k line trjatomnum trjcollection atomid_conv +clear BOinform lineofmolecule numofmolecule tarBOrow tarmatchcol tartrjdata trjrow trjdataname trjlength trjmax trjmin +clear trjmod trjnnum trjnum trjones trjper trjrawdata trjreadline trjstep MOLE lib base eleswapans eleswap boxsize formatout + \ No newline at end of file