atat-mirror/src/nnshell.c++

87 lines
2.8 KiB
C++

#include <fstream.h>
#include "xtalutil.h"
#include "getvalue.h"
#include "parse.h"
#include "version.h"
void list_atom(ostream &file, const Array<Arrayint> &labellookup, const Array<AutoString> &label, int thetype) {
for (int i=0; i<labellookup(thetype).get_size(); i++) {
if (i>0) {cout << ",";}
cout << label(labellookup(thetype)(i));
}
}
int main(int argc, char *argv[]) {
char *latfilename="lat.in";
int listbond=0;
Real user_rad=0.0;
int sigdig=5;
int dummy;
AskStruct options[]={
{"","Nearest Neighbor Shell " MAPS_VERSION ", by Axel van de Walle and Paul Dalach",TITLEVAL,NULL},
{"-l","Input file defining the lattice (default: lat.in)",STRINGVAL,&latfilename},
{"-lb","List bond types in given structures",BOOLVAL,&listbond},
{"-sig","Number of significant digits to print in output files",INTVAL,&sigdig},
{"-r","Define a max radius for neighbor search",REALVAL,&user_rad},
{"-d","Use all default values",BOOLVAL,&dummy}
};
if (!get_values(argc,argv,countof(options),options)) {
display_help(countof(options),options);
return 1;
}
Structure lat;
Array<Arrayint> labellookup;
Array<AutoString> label;
rMatrix3d axes;
if (strcmp(latfilename,"-")==0) {
parse_lattice_file(&lat.cell, &lat.atom_pos, &lat.atom_type, &labellookup, &label, cin, &axes);
}
else {
ifstream latfile(latfilename);
if (!latfile) ERRORQUIT("Unable to open lattice file");
parse_lattice_file(&lat.cell, &lat.atom_pos, &lat.atom_type, &labellookup, &label, latfile, &axes);
}
wrap_inside_cell(&lat.atom_pos,lat.atom_pos,lat.cell);
cout.setf(ios::fixed);
cout.precision(sigdig);
if (listbond) {
for (int at=0; at<lat.atom_type.get_size(); at++) {
Real rnn=find_1nn_radius(lat,at);
AtomPairIterator p(lat.cell,lat.atom_pos(at),lat.atom_pos);
for ( ; p.length()<=rnn+zero_tolerance; p++) {
for (int i=0; i<2; i++) {
if (i>0) {cout <<"-";}
int thetype=lat.atom_type(which_atom(lat.atom_pos,p(i),!lat.cell));
list_atom(cout,labellookup,label,thetype);
}
cout << endl;
}
}
}
else {
for (int at=0; at<lat.atom_type.get_size(); at++) {
int nbnn=0;
cout << (!axes)*lat.atom_pos(at) << endl;
Real rnn;
if(user_rad>0) {rnn=user_rad;}
else {rnn=find_1nn_radius(lat,at);}
cout << rnn << endl;
list_atom(cout,labellookup,label,lat.atom_type(at));
cout << endl;
AtomPairIterator p(lat.cell,lat.atom_pos(at),lat.atom_pos);
while (p.length()<=rnn+zero_tolerance) {
cout << " " << (!axes)*(p(1)-p(0)) << " " << norm(p(1)-p(0)) << " ";
int thetype=lat.atom_type(which_atom(lat.atom_pos,p(1),!lat.cell));
list_atom(cout,labellookup,label,thetype);
cout << endl;
nbnn++;
p++;
}
cout << " " << nbnn << endl << endl;
}
}
}