Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add -csa flag for subsampling pore surfaces #3

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

add -csa flag for subsampling pore surfaces #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
138 changes: 133 additions & 5 deletions area_and_volume.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,133 @@ void reportPointsVisIT(ostream &output, vector<Point> axsPoints, vector<int> axs
output << coords[0] << " " << coords[1] << " " << coords[2] << " 0 n " << inaxsPPIDs[i] << "\n";
}
}
void reportPointsCIF(ostream &output, vector<Point> axsPoints, vector<int> axsPChIDs, vector<Point> inaxsPoints, vector<int> inaxsPPIDs, ATOM_NETWORK *cell,double totASA,double totNASA){


string formula = "pore-vis";
output << "data_" << formula << endl;
output << "#******************************************" << endl;
output << "#" << endl;
output << "# CIF file created by Zeo++" << endl;
output << "# Zeo++ is an open source package to" << endl;
output << "# analyze microporous materials" << endl;
output << "#" << endl;
output << "#*******************************************" << "\n\n";
output << "_cell_length_a\t\t" << cell->a << " " << endl; // removed (0)
output << "_cell_length_b\t\t" << cell->b << " " << endl;
output << "_cell_length_c\t\t" << cell->c << " " << endl;
output << "_cell_angle_alpha\t\t" << cell->alpha << " " << endl;
output << "_cell_angle_beta\t\t" << cell->beta << " " << endl;
output << "_cell_angle_gamma\t\t" << cell->gamma << " \n\n";
output << "_symmetry_space_group_name_H-M\t\t" << "'P1'" << endl;
output << "_symmetry_Int_Tables_number\t\t" << "1" << endl;
output << "_symmetry_cell_setting\t\t";

//Determine the Crystal System
if (cell->alpha == 90 && cell->beta == 90 && cell->gamma == 90){
if (cell->a == cell->b || cell->b == cell->c || cell->a == cell->c){
if (cell->a == cell->b && cell->b == cell->c){
output << "Isometric\n" << endl;
}
else {
output << "Tetragonal\n" << endl;
}
}
else{
output << "Orthorhombic\n" << endl;
}
}
else if(cell->alpha == cell->beta || cell->beta == cell->gamma || cell->alpha == cell->gamma){
output << "Monoclinic\n" << endl;
}
else{
output << "Triclinic\n" << endl;
}

output << "loop_" << endl;
output << "_symmetry_equiv_pos_as_xyz" << endl;
output << "'+x,+y,+z'\n" << endl;
output << "loop_" << endl;
output << "_atom_site_label" << endl;
output << "_atom_site_type_symbol" << endl;
output << "_atom_site_fract_x" << endl;
output << "_atom_site_fract_y" << endl;
output << "_atom_site_fract_z" << endl;
output << "_atom_site_calc_flag" << endl; // channel id
output << "_atom_site_description" << endl; // accessible:0, inaccessible 1

// Performing subsampling of points //
vector<Point> sub_axsPoints = vector<Point> (); // stores accessible points
vector<Point> sub_inaxsPoints = vector<Point> (); // stores accessible pointsA
int numAxsSampled=(int)(0.3*totASA);
int numInAxsSampled=(int)(0.3*totNASA);
cout << "Subsampling points for the cif output.." <<"\n";
if (axsPoints.size()!=0){
sub_axsPoints.push_back(axsPoints.at(0)); // adding one point to the set
double rcut=1.3;
int count_fail=0;
while (sub_axsPoints.size()< numAxsSampled){
int selected = rand() % axsPoints.size();
Point coords_selected = axsPoints.at(selected);
double min_d=10000.0;
for(unsigned int i = 0; i < sub_axsPoints.size(); i++)
{
Point coords = sub_axsPoints.at(i);
double d_pair = cell->calcDistanceABC( coords_selected[0] , coords_selected[1] , coords_selected[2] , coords[0],coords[1],coords[2]);
if (d_pair < min_d) min_d=d_pair;
}
if (min_d > rcut){
sub_axsPoints.push_back(coords_selected);
} else{
if (count_fail == 2000){
count_fail=0;
rcut-=0.1;
if (rcut <= 0.7) break;
}else {
count_fail+=1;
}
}
}
}
if (inaxsPoints.size()!=0){
sub_inaxsPoints.push_back(inaxsPoints.at(0));
double rcut=1.3;
int count_fail=0;
while (sub_inaxsPoints.size()< numInAxsSampled){
int selected = rand() % inaxsPoints.size();
Point coords_selected = inaxsPoints.at(selected);
double min_d=10000.0;
for(unsigned int i = 0; i < sub_inaxsPoints.size(); i++)
{
Point coords = sub_inaxsPoints.at(i);
double d_pair = cell->calcDistanceABC( coords_selected[0] , coords_selected[1] , coords_selected[2] , coords[0],coords[1],coords[2]);
if (d_pair < min_d) min_d=d_pair;
}
if (min_d > rcut){
sub_inaxsPoints.push_back(coords_selected);
} else{
if (count_fail == 2000){
count_fail=0;
rcut-=0.1;
if (rcut <= 0.7) break;
}else {
count_fail+=1;
}
}
}
}


for(unsigned int i = 0; i < sub_axsPoints.size(); i++){
Point coords = sub_axsPoints.at(i);
output << "A " << " A " << coords[0] << " " << coords[1] << " " << coords[2] << " "<< axsPChIDs[i] << " 0 " <<"\n";
}
for(unsigned int i = 0; i < sub_inaxsPoints.size(); i++){
Point coords = sub_inaxsPoints.at(i);
output << "N "<< " N " << coords[0] << " " << coords[1] << " " << coords[2] << " "<< inaxsPPIDs[i] <<" 1 " <<"\n";
}

}


/* To adjust the sampling point to minimize its distance with the central atom, shift the point by the reverse of the
Expand Down Expand Up @@ -1674,7 +1801,7 @@ void getLocalSubstructures(char* name, double probeRad, double local_substructur
}


double calcASA(ATOM_NETWORK *hiaccatmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe_chan, double r_probe, double rho_crystal, int numSamples, bool excludePockets, ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool ExtendedOutputFlag){
double calcASA(ATOM_NETWORK *hiaccatmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe_chan, double r_probe, double rho_crystal, int numSamples, bool excludePockets, ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool CIFOFlag, bool ExtendedOutputFlag){

ATOM_NETWORK *atmnet; // pointed to the analyzed (original atomsnet)

Expand Down Expand Up @@ -1835,7 +1962,8 @@ double calcASA(ATOM_NETWORK *hiaccatmnet, ATOM_NETWORK *orgatmnet, bool highAccu
LiverpoolInaxsPoints.push_back(atmnet->xyz_to_abc(inaxsPoints[j]));
};
// reportPointsVisIT(output, LiverpoolAxsPoints, LiverpoolInaxsPoints);
reportPointsVisIT(output, LiverpoolAxsPoints, axsPointsChannelIDs, LiverpoolInaxsPoints, inaxsPointsPocketIDs);
//reportPointsVisIT(output, LiverpoolAxsPoints, axsPointsChannelIDs, LiverpoolInaxsPoints, inaxsPointsPocketIDs);
reportPointsCIF(output, LiverpoolAxsPoints, axsPointsChannelIDs, LiverpoolInaxsPoints, inaxsPointsPocketIDs, atmnet,totalSA, totalSA_inaxs);
};
};
//reportResampledPoints(output, resampledInfo);
Expand Down Expand Up @@ -2455,8 +2583,8 @@ bool inside = overlaps; // Only want to visualize points that are not inside sam
*/


double calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe, double rho_crystal, int numSamples, bool excludePockets, ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool ExtendedOutputFlag){
return calcASA(atmnet, orgatmnet, highAccuracy, r_probe, r_probe, rho_crystal, numSamples, excludePockets, output, filename, visualize, VisITflag, LiverpoolFlag, ExtendedOutputFlag);
double calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe, double rho_crystal, int numSamples, bool excludePockets, ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool CIFOFlag, bool ExtendedOutputFlag){
return calcASA(atmnet, orgatmnet, highAccuracy, r_probe, r_probe, rho_crystal, numSamples, excludePockets, output, filename, visualize, VisITflag, LiverpoolFlag,CIFOFlag, ExtendedOutputFlag);
}

/// Cython wrapper to calcASA where visualization flags are ignored
Expand All @@ -2465,7 +2593,7 @@ string calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy,
stringstream output;
string filename = "No filename";
double rho_crystal = calcDensity(atmnet);
double sa = calcASA(atmnet, orgatmnet, highAccuracy, r_probe_chan, r_probe, rho_crystal, numSamples, excludePockets, output, (char *)filename.data(), false, false, false, ExtendedOutputFlag);
double sa = calcASA(atmnet, orgatmnet, highAccuracy, r_probe_chan, r_probe, rho_crystal, numSamples, excludePockets, output, (char *)filename.data(), false, false, false,false, ExtendedOutputFlag);
return output.str();
}

Expand Down
5 changes: 3 additions & 2 deletions area_and_volume.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ double calcDensity(ATOM_NETWORK *atmnet);
* inaccessible points are colored red.*/
void reportPoints(std::ostream &output, std::vector<Point> axsPoints, std::vector<Point> inaxsPoints);
void reportPointsVisIT(std::ostream &output, std::vector<Point> axsPoints, std::vector<Point> inaxsPoints);
void reportPointsCIF(ostream &output, vector<Point> axsPoints, vector<int> axsPChIDs, vector<Point> inaxsPoints, vector<int> inaxsPPIDs,ATOM_NETWORK *cell,double totSA, double totNASA);

/* Print the coordinates contained in the provided vectors in a manner that they can be
* displayed using ZeoVis and its VMD interface. Useful for debugging errors in the
Expand Down Expand Up @@ -57,8 +58,8 @@ double calcAV(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy,
// For python interface
std::string calcAV(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe_chan, double r_probe, int numSamples, bool excludePockets, double low_dist_cutoff, double high_dist_cutoff);

double calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe_chan, double r_probe, double rho_crystal, int numSamples, bool excludePockets, std::ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool ExtendedOutputFlag);
double calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe, double rho_crystal, int numSamples, bool excludePockets, std::ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool ExtendedOutputFlag);
double calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe_chan, double r_probe, double rho_crystal, int numSamples, bool excludePockets, std::ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool CIFOFlag ,bool ExtendedOutputFlag);
double calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe, double rho_crystal, int numSamples, bool excludePockets, std::ostream &output, char *filename, bool visualize, bool VisITflag, bool LiverpoolFlag, bool CIFOFlag,bool ExtendedOutputFlag);
// For python interface
std::string calcASA(ATOM_NETWORK *atmnet, ATOM_NETWORK *orgatmnet, bool highAccuracy, double r_probe_chan, double r_probe, int numSamples, bool excludePockets, bool ExtendedOutputFlag);

Expand Down
19 changes: 12 additions & 7 deletions main.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -521,7 +521,7 @@ int main(int argc, char * argv[]){
double volume = calcDeterminant(atmnet.ucVectors); // Unit cell volume/Units of A^3
int numSamples = (int)(volume*numSamplesAV);
calcAV(&atmnet, &orgAtomnet, highAccuracy, probe_radius, probe_radius, numSamples, true, output, (char *)filename.c_str(), 0, 0, 0, 0, -1,-1, 0);
calcASA(&atmnet, &orgAtomnet, highAccuracy, probe_radius, probe_radius, calcDensity(&atmnet), numSamplesASA, true, output, (char *)filename.c_str(), 0, 0, 0, 0);
calcASA(&atmnet, &orgAtomnet, highAccuracy, probe_radius, probe_radius, calcDensity(&atmnet), numSamplesASA, true, output, (char *)filename.c_str(), 0, 0, 0, 0,0);

for(unsigned int i = 0; i < pores.size(); i++)
if(pores[i].dimensionality>0) pores[i].printPoreSummary(output, &atmnet); // Channels
Expand Down Expand Up @@ -550,16 +550,21 @@ int main(int argc, char * argv[]){
// calculateAverageGrid(&atmnet, filename_InputData, filename_Gaussian_cube, angstrom_to_bohr, useMass);
analyzePoreInfoFiles(&atmnet, filename_InputData, filename_output);
}
else if((command[0].compare("-sa") == 0) || (command[0].compare("-saex") == 0) || (command[0].compare("-zsa") == 0) || (command[0].compare("-vsa") == 0) || (command[0].compare("-lsa")==0)){
bool visualize = (command[0].compare("-zsa") == 0 || command[0].compare("-vsa") == 0 || command[0].compare("-lsa")==0);
bool visVisITflag = (command[0].compare("-vsa")==0 || command[0].compare("-lsa")==0);
bool LiverpoolFlag = (command[0].compare("-lsa")==0);
else if((command[0].compare("-sa") == 0) || (command[0].compare("-saex") == 0) || (command[0].compare("-zsa") == 0) || (command[0].compare("-vsa") == 0) || (command[0].compare("-lsa")==0)||(command[0].compare("-csa")==0)){
bool visualize = (command[0].compare("-zsa") == 0 || command[0].compare("-vsa") == 0 || command[0].compare("-lsa")==0||(command[0].compare("-csa")==0));
bool visVisITflag = (command[0].compare("-vsa")==0 || command[0].compare("-lsa")==0||(command[0].compare("-csa")==0));
bool LiverpoolFlag = (command[0].compare("-lsa")==0||(command[0].compare("-csa")==0));
bool CIFOFlag = ((command[0].compare("-csa")==0));
bool ExtendedOutputFlag = (command[0].compare("-saex") == 0);
string suffix;
if(visualize)
{
if(visVisITflag) {
if(LiverpoolFlag) suffix = ".lsa"; else suffix = ".vsa";
if(LiverpoolFlag && CIFOFlag)
suffix = "_SA.cif";
else if (LiverpoolFlag)
suffix = ".lsa";
else suffix = ".vsa";
} else {suffix = ".zsa";};
}
else
Expand All @@ -572,7 +577,7 @@ int main(int argc, char * argv[]){

fstream output;
output.open(filename.data(), fstream::out);
calcASA(&atmnet, &orgAtomnet, highAccuracy, chan_radius, probe_radius, calcDensity(&atmnet), numSamples, true, output, (char *)filename.data(), visualize, visVisITflag, LiverpoolFlag, ExtendedOutputFlag);
calcASA(&atmnet, &orgAtomnet, highAccuracy, chan_radius, probe_radius, calcDensity(&atmnet), numSamples, true, output, (char *)filename.data(), visualize, visVisITflag, LiverpoolFlag,CIFOFlag, ExtendedOutputFlag);
output.close();
}
else if( (command[0].compare("-vol") == 0) || (command[0].compare("-zvol") == 0) || (command[0].compare("-vvol") == 0) || (command[0].compare("-lvol") == 0) || (command[0].compare("-volpo") == 0)) {
Expand Down
2 changes: 1 addition & 1 deletion zeojobs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -527,7 +527,7 @@ class zeoJob {
double volume = calcDeterminant(Material.atmnet.ucVectors); // Unit cell volume/Units of A^3
int numSamples = (int)(volume*numSamplesAV);
calcAV(&(Material.atmnet), &(Material.orgAtomnet), Material.highAccuracy, probeRadius, probeRadius, numSamples, true, output, (char *)filename.c_str(), 0, 0, 0, 0, -1,-1, 0);
calcASA(&(Material.atmnet), &(Material.orgAtomnet), Material.highAccuracy, probeRadius, probeRadius, calcDensity(&(Material.atmnet)), numSamplesASA, true, output, (char *)filename.c_str(), 0, 0, 0, 0);
calcASA(&(Material.atmnet), &(Material.orgAtomnet), Material.highAccuracy, probeRadius, probeRadius, calcDensity(&(Material.atmnet)), numSamplesASA, true, output, (char *)filename.c_str(), 0, 0, 0, 0,0);

for(unsigned int i = 0; i < pores.size(); i++)
if(pores[i].dimensionality>0) pores[i].printPoreSummary(output, &(Material.atmnet)); // Channels
Expand Down