diff --git a/C12E2_work/#C12E2FullAnalysis.cxx# b/C12E2_work/#C12E2FullAnalysis.cxx# new file mode 100755 index 00000000..f3915efb --- /dev/null +++ b/C12E2_work/#C12E2FullAnalysis.cxx# @@ -0,0 +1,999 @@ +/* + + ------------------------------------------------------------------------ + | Analysis of C12E2 Mixed Bilayers | + ------------------------------------------------------------------------ + + Author = "Sang Young Noh" + + Version = "0.0.1" + + Updated = 11/09/2019 + + The papers referenced for this work is: + + 1. "F. C. MacKintosh and S. A. Safran, Phase Separation and curvature of + bilayer membranes, Physical Review E, 47, 2, 1993" + + 2. "R. Lipowsky, Domain-induced budding of fluid membranes, Biophysical + Society, 64, 1993, 1133 - 1138" + + 3. "J. Wolff and S. Komura and D. Andelman, Budding of domains in mixed bilayer + domains, 91, Physical Review E, 91, 012708, 2015" + + 4. S. Katira and K. K. Madadapu and S. Vaikuntanathan and B. Smit and D. + Chandler, Pre-transition effects mediate forces of assembly between + transmembrane proteins, elife, 2016, 5, e13150 + + 5. Physical Considerations of the Organization of Inclusions in Lipid Bilayer + Systems, Spring 2015, Shachi Katira + +*/ + +#include +#include +#include +#include +#include +#include // Library for complex numbers +#include +#include // std::multiplies +#include +#include +#include +#include // std::adjacent_differenc#include +#include +#include +#include + +#define _USE_MATH_DEFINES + +// Constants for use later with analysis +const int NumberOfPolymers = + 998; // The number of polymers of each type - C12E2 or mimic +const int numberofatoms = 71313; // Total number of beads in the simulation +const int indexCG = 7; +int numberofSS = 100; /*The number of screenshots in the dump file*/ +const int boxdim = 3; + +typedef struct { + int index[7]; +} C12E2_skeleton; + +typedef struct { + int a; + int b; + double x; + double y; + double z; +} inputCoord; + +typedef struct { // Used to input the center of masses for each lipid + int index; + double x; + double y; + double z; +} COMstruct; + +typedef struct { + double dist; + int topPhiC12E2Count; + int topPhiC12E2MCount; + int botPhiC12E2Count; + int botPhiC12E2MCount; +} phiStruct; + +typedef struct { + double phip; + double phim; + std::vector phipVec; + std::vector phimVec; +} phipm; + +typedef struct { + double X, Y, Z; + double orderphobicVal; + std::vector surrounding_residue_indices; +} OPVal; + +typedef struct { // Used to identify the group and distance to compute the + // orderphobic effect + int index; + double dist; + double selfXcoord; + double selfYcoord; + double selfZcoord; + double Xcoord; + double Ycoord; + double Zcoord; + +} OPHstruct; + +double calcAngle(std::vector *OPHInput) { + + double refVector[2]; + double dot; + double newx; + double newy; + double u, v; + + std::complex angle; + + double phi; + + refVector[0] = 10.0; + + refVector[1] = 10.0; + + std::complex mycomplex(0.0, 0.0); + std::complex sixth((double)1 / 6, 0.0); + + double absval; + + for (unsigned int i = 0; i < OPHInput->size(); ++i) { + + newx = OPHInput->at(i).Xcoord; + + newy = OPHInput->at(i).Ycoord; + + u = pow((pow(refVector[0], 2) + pow(refVector[1], 2)), 0.5); + + v = pow((pow(newx, 2) + pow(newy, 2)), 0.5); + + dot = (newx * refVector[0]) + (newy * refVector[1]); + + std::complex placeholder( + 0, (6 * 180 * M_PI * acos((dot / (u * v))))); // TODO + + angle = exp(placeholder); + mycomplex += angle; + + // std::cout << i << " " << absval << std::endl; + } + mycomplex = sixth * mycomplex; + absval = pow(abs(mycomplex), 2); + // std::cout << absval << std::endl; + return absval; +} + +double trueDist(double *COM1x, double *COM1y, double *COM1z, double *COM2x, + double *COM2y, double *COM2z) { + double dist = pow((pow(*COM1x - *COM2x, 2) + pow(*COM1y - *COM2y, 2) + + pow(*COM1z - *COM2z, 2)), + 0.5); + return dist; +} + +// Function used for sorting vector of structs +bool compareByIndex(const inputCoord &a, const inputCoord &b) { + return a.a < b.a; +} + +bool compareByIndexOPh(const OPHstruct &a, const OPHstruct &b) { + return a.dist < b.dist; +} + +void CenterOfMass(std::vector *vec1, std::vector *vec2, + std::vector> *inputVec, + std::vector *COM1, std::vector *COM2, + std::vector> *C12E2final, + std::vector> *C12E2Mfinal) { + // CenterOfMass(&C12E2IndexVector, &C12E2MIndexVector, &inputTotal, + // &C12E2COM, &C12E2MCOM, &C12E2TotalCOMArray, &C12E2MTotalCOMArray); + + double C12E2comX, C12E2McomX; // X coordinate COM + double C12E2comY, C12E2McomY; // Y coordinate COM + double C12E2comZ, C12E2McomZ; // Z coordinate cCOM + + COMstruct C12E2input; // struct to store COM coordinates for C12E2 + COMstruct C12E2Minput; // struct to store cOM coordinates for C12E2M + + for (unsigned int i = 0; i <= inputVec->size() - 1; ++i) { + + for (unsigned int index = 0; index <= vec1->size() - 1; ++index) { + + C12E2comX = (inputVec->at(i).at(vec1->at(index)).x + + inputVec->at(i).at((vec1->at(index)) + 1).x + + inputVec->at(i).at((vec1->at(index)) + 2).x + + inputVec->at(i).at((vec1->at(index)) + 3).x + + inputVec->at(i).at((vec1->at(index)) + 4).x + + inputVec->at(i).at((vec1->at(index)) + 5).x + + inputVec->at(i).at((vec1->at(index)) + 6).x) / + 7.0; + + C12E2comY = (inputVec->at(i).at(vec1->at(index)).y + + inputVec->at(i).at((vec1->at(index)) + 1).x + + inputVec->at(i).at((vec1->at(index)) + 2).y + + inputVec->at(i).at((vec1->at(index)) + 3).y + + inputVec->at(i).at((vec1->at(index)) + 4).y + + inputVec->at(i).at((vec1->at(index)) + 5).y + + inputVec->at(i).at((vec1->at(index)) + 6).y) / + 7.0; + + C12E2comZ = (inputVec->at(i).at(vec1->at(index)).z + + inputVec->at(i).at((vec1->at(index)) + 1).z + + inputVec->at(i).at((vec1->at(index)) + 2).z + + inputVec->at(i).at((vec1->at(index)) + 3).z + + inputVec->at(i).at((vec1->at(index)) + 4).z + + inputVec->at(i).at((vec1->at(index)) + 5).z + + inputVec->at(i).at((vec1->at(index)) + 6).z) / + 7.0; + + C12E2input.index = vec1->at(index); + C12E2input.x = C12E2comX; + C12E2input.y = C12E2comY; + C12E2input.z = C12E2comZ; + COM1->push_back(C12E2input); + } + + for (unsigned int index = 0; index <= vec2->size() - 1; ++index) { + + C12E2McomX = (inputVec->at(i).at(vec2->at(index)).x + + inputVec->at(i).at((vec2->at(index)) + 1).x + + inputVec->at(i).at((vec2->at(index)) + 2).x + + inputVec->at(i).at((vec2->at(index)) + 3).x + + inputVec->at(i).at((vec2->at(index)) + 4).x + + inputVec->at(i).at((vec2->at(index)) + 5).x + + inputVec->at(i).at((vec2->at(index)) + 6).x) / + 7.0; + + C12E2McomY = (inputVec->at(i).at(vec2->at(index)).y + + inputVec->at(i).at((vec2->at(index)) + 1).x + + inputVec->at(i).at((vec2->at(index)) + 2).y + + inputVec->at(i).at((vec2->at(index)) + 3).y + + inputVec->at(i).at((vec2->at(index)) + 4).y + + inputVec->at(i).at((vec2->at(index)) + 5).y + + inputVec->at(i).at((vec2->at(index)) + 6).y) / + 7.0; + + C12E2McomZ = (inputVec->at(i).at(vec2->at(index)).z + + inputVec->at(i).at((vec2->at(index)) + 1).z + + inputVec->at(i).at((vec2->at(index)) + 2).z + + inputVec->at(i).at((vec2->at(index)) + 3).z + + inputVec->at(i).at((vec2->at(index)) + 4).z + + inputVec->at(i).at((vec2->at(index)) + 5).z + + inputVec->at(i).at((vec2->at(index)) + 6).z) / + 7.0; + + C12E2Minput.index = vec2->at(index); + C12E2Minput.x = C12E2comX; + C12E2Minput.y = C12E2comY; + C12E2Minput.z = C12E2comZ; + // TODO + COM1->push_back(C12E2Minput); + } + + C12E2final->push_back(*COM1); + // C12E2final->push_back(*COM2); + // TODO + C12E2Mfinal->push_back(*COM2); + COM1->clear(); + COM2->clear(); + } +} + +bool sortbysec_int(const std::pair &a, const std::pair &b) { + return (a.second < b.second); +} + +bool sortbysec_double(const std::pair &a, + const std::pair &b) { + return (a.second < b.second); +} + +class compute { +public: + compute(){}; // Default constructor + compute(std::string file){}; // constructor reading file + + void storeFile() { + + boost::progress_display show_progress(numberofSS); + + // open file for reading + ipf = fopen("dump.mixed2", "r"); // Needs correction + // check for error opening file */ + if (ipf == NULL) { + std::cout << "Error opening file\n"; + exit(1); + } + // loop over the values + for (int SSno = 0; SSno < numberofSS; SSno++) { + + l = 0; + n = 0; + index = 0; + // Make sure to clear all vectors before going on to the next snapshot + inputVector.clear(); + for (int k = 0; k < nlines; k++) { + // get a line from the file + // fgets() returns NULL when it reaches an error or end of file + + fgets(line, sizeof(line), ipf); + // while (fgets(line,sizeof(line),ipf) != NULL) { + + if (l < 5) { + /* We are doing nothing */ + // printf("%s l= %d", line,l); + l++; + } else if ((l > 4 && l < 8)) { + /*We are scanning the bit with just the box parameters*/ + sscanf(line, "%lf %lf", &box1, &box2); + // printf("%lf %lf \n",box1,box2 ); + boxlength[l - 5] = box2 - box1; + l++; + } else if (l == 8) { + // printf(" **** l= %d \n",l); + /* We are doing nothing */ + l++; + } else { + // printf(" ***** l = %d \n",l ); + /* convert the text to numbers */ + + sscanf(line, "%d %d %lf %lf %lf", &index, &atomtype, &x, &y, &z); + // std::cout << index << " " << atomtype << " " << x << " " << y << " + // " << z; a.push_back(std::make_pair(index, index)); // Push back + // indices b.push_back(std::make_pair(index, atomtype)); // Push back + // atomtypes xco.push_back(std::make_pair(index, x*boxlength[0])); // + // Push back boxlengths - x coordinates + // yco.push_back(std::make_pair(index, y*boxlength[1])); // Push back + // boxlengths - y coordinates zco.push_back(std::make_pair(index, + // z*boxlength[2])); // Push back boxlengths - z coordinates + inputline.a = index; + inputline.b = atomtype; + inputline.x = x * boxlength[0]; + inputline.y = y * boxlength[1]; + inputline.z = z * boxlength[2]; + inputVector.push_back(inputline); + n++; + l++; + // printf("%d %d %lf %lf %lf\n",index,atomtype,x,y,z); + } + } + inputTotal.push_back(inputVector); + ++show_progress; + } + } + + void sortVectors() { // Sort vector values based on indices of the dump files + for (unsigned int i = 0; i < inputTotal.size(); ++i) { + sort(inputTotal[i].begin(), inputTotal[i].end(), compareByIndex); + } + } + + void printVectorElements() { + for (unsigned int i = 0; i < inputTotal.size(); ++i) { + for (unsigned int j = 0; j < inputTotal[1].size(); ++j) { + // std::cout << i << " " << j << " " << " " << inputTotal[i][j].a << " " + // << inputTotal[i][j].b << " " << inputTotal[i][j].x << " " << + // inputTotal[i][j].y << " " << inputTotal[i][j].z << std::endl; + } + } + } + + /* + void CompositionProfile() { + //Calculating the composition profile + for (int i = 0; i <= sizeof(C12E2_struct)/sizeof(C12E2_struct[1])-1; i++) { + // TODO + + //double truedist(double COM1x, double COM1y, double COM1z, double COM2x, + double COM2y, double COM2z) { + // See if the distance between the nanopparticle is smaller than 15 + //std::cout << zco[71312] << std::endl; + + if (trueDist(xco[71312], yco[71312], zco[71312], headGroupC12_E2xCOM[i], + headGroupC12_E2yCOM[i], headGroupC12_E2zCOM[i]) <= 25.000) { + // Check if the c12e2 molecule is on the top layer or the bottom layer + if (headGroupC12_E2zCOM[i] > zco[71312]) { // If on the top layer + tophead++; + } + else if (headGroupC12_E2zCOM[i] < zco[71312]) { // If on the bottom + layer downhead++; + } + } + + // This time, working with mimics + if(trueDist(xco[71312], yco[71312], zco[71312], headGroupMIMCxCOM[i], + headGroupMIMICyCOM[i], headGroupMIMICzCOM[i] ) <= 25.000) { if + (headGroupMIMICzCOM[i] > zco[71312]) { mimictophead++; + } + else if (headGroupMIMICzCOM[i] < zco[71312]) { + mimicdownhead++; + } + } + } + + + //std::cout << SSno << " " << tophead/(tophead + mimictophead) << " " << + mimictophead/(tophead + mimictophead) << " " << downhead/(downhead + + mimicdownhead) << " " << mimicdownhead/(downhead + mimicdownhead) << " " << + (tophead-mimictophead)/(tophead + mimictophead) << " " << + (downhead-mimicdownhead)/(downhead + mimicdownhead) << " " << + ((tophead-mimictophead)/(tophead + mimictophead) + + (downhead-mimicdownhead)/(downhead + mimicdownhead))/2- << " " << + ((tophead-mimictophead)/(tophead + mimictophead) - + (downhead-mimicdownhead)/(downhead + mimicdownhead))/2 << std::endl; + } + */ + + void check() { + for (unsigned int i = 0; i < inputTotal.size(); ++i) { + for (unsigned int j = 0; j <= inputTotal[1].size(); j++) { + if (inputTotal[i][j].b == 7) { + std::cout << inputTotal[i][j].a << " " << inputTotal[i][j].b + << std::endl; + std::cout << inputTotal[i][j + 1].a << " " << inputTotal[i][j + 1].b + << std::endl; + std::cout << inputTotal[i][j + 2].a << " " << inputTotal[i][j + 2].b + << std::endl; + std::cout << inputTotal[i][j + 3].a << " " << inputTotal[i][j + 3].b + << std::endl; + std::cout << inputTotal[i][j + 4].a << " " << inputTotal[i][j + 4].b + << std::endl; + std::cout << inputTotal[i][j + 5].a << " " << inputTotal[i][j + 5].b + << std::endl; + std::cout << inputTotal[i][j + 6].a << " " << inputTotal[i][j + 6].b + << std::endl; + } else if (inputTotal[i][j].b == 13) { + std::cout << inputTotal[i][j].a << " " << inputTotal[i][j].b + << std::endl; + std::cout << inputTotal[i][j + 1].a << " " << inputTotal[i][j + 1].b + << std::endl; + std::cout << inputTotal[i][j + 2].a << " " << inputTotal[i][j + 2].b + << std::endl; + std::cout << inputTotal[i][j + 3].a << " " << inputTotal[i][j + 3].b + << std::endl; + std::cout << inputTotal[i][j + 4].a << " " << inputTotal[i][j + 4].b + << std::endl; + std::cout << inputTotal[i][j + 5].a << " " << inputTotal[i][j + 5].b + << std::endl; + std::cout << inputTotal[i][j + 6].a << " " << inputTotal[i][j + 6].b + << std::endl; + } + } + } + } + + void headGroupVectorFormation() { + for (unsigned int j = 0; j <= inputTotal[1].size(); j++) { + if (inputTotal[1][j].b == 7) { + C12E2IndexVector.push_back( + inputTotal[1][j].a); // push back C12E2 bead 7 indices (headgroups) + } else if (inputTotal[1][j].b == 13) { + C12E2MIndexVector.push_back( + inputTotal[1][j].a); // push back C12E2 bead 7 indices (headgroups) + } + } + + for (unsigned int index = 0; index <= C12E2IndexVector.size(); ++index) { + std::cout << C12E2IndexVector[index] << std::endl; + } + + for (unsigned int index = 0; index <= C12E2MIndexVector.size(); ++index) { + std::cout << C12E2MIndexVector[index] << std::endl; + } + } + + void ComputePhi() { // Computes the phi, or the mismatch between the bilayer + // leaflets around the NP + // std::cout << C12E2IndexVector.size() << " " << C12E2MIndexVector.size() + // << " " << inputTotal.size() << std::endl; + + /* + Following Python's comment system """ """ + + We are explcitily ignoring flip-flops in the case of this study + */ + + for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { + if (inputTotal[1][C12E2IndexVector[index]].z > 120.0) { + + topC12E2Index.push_back(C12E2IndexVector[index]); + } else if (inputTotal[1][C12E2IndexVector[index]].z < 120.0) { + botC12E2Index.push_back(C12E2IndexVector[index]); + } + } + + for (unsigned int index = 0; index < C12E2MIndexVector.size(); ++index) { + if (inputTotal[1][C12E2MIndexVector[index]].z > 120.0) { + topC12E2MIndex.push_back(C12E2MIndexVector[index]); + } else if (inputTotal[1][C12E2MIndexVector[index]].z < 120.0) { + botC12E2MIndex.push_back(C12E2MIndexVector[index]); + } + } + + std::cout << "check" << std::endl; + + for (std::vector::iterator it = topC12E2Index.begin(); + it != topC12E2Index.end(); it++) { + // std::cout << *it << " " << std::endl; + // std::cout << topC12E2Index.size() << std::endl; + } + + for (std::vector::iterator it = botC12E2Index.begin(); + it != botC12E2Index.end(); it++) { + // std::cout << *it << " " << std::endl; + // std::cout << botC12E2Index.size() << std::endl; + } + + std::cout << topC12E2Index.size() << std::endl; + std::cout << botC12E2Index.size() << std::endl; + + std::cout << topC12E2MIndex.size() << std::endl; + std::cout << botC12E2MIndex.size() << std::endl; + phiStruct phitemplate; + std::cout << "check" << std::endl; + + for (unsigned int i = 0; i <= inputTotal.size() - 1; ++i) { + + for (unsigned phiIndex = 0; phiIndex <= 100; ++phiIndex) { + phitemplate.dist = phiIndex; + phitemplate.topPhiC12E2Count = 0; + phitemplate.botPhiC12E2Count = 0; + phitemplate.topPhiC12E2MCount = 0; + phitemplate.botPhiC12E2MCount = 0; + phiCount.push_back(phitemplate); + } + + NPX = inputTotal[i][71312].x; // x coordinate of the NP + NPY = inputTotal[i][71312].y; // y coordinate of the NP + NPZ = inputTotal[i][71312].z; // z coordinate of the NP + + int topC12E2Counter, botC12E2Counter, topC12E2MCounter, botC12E2MCounter; + + for (unsigned int dist = 0; dist <= 100; ++dist) { + + for (unsigned int index = 0; index < topC12E2Index.size(); index++) { + DistVec1 = + trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][topC12E2Index[index]].x, + &inputTotal[i][topC12E2Index[index]].y, + &inputTotal[i][topC12E2Index[index]].z); + if (DistVec1 <= dist + 0.5 && DistVec1 >= dist - 0.5) { + phiCount[dist].topPhiC12E2Count += 1; + } + } + + for (unsigned int index = 0; index < botC12E2Index.size(); index++) { + DistVec2 = + trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][botC12E2Index[index]].x, + &inputTotal[i][botC12E2Index[index]].y, + &inputTotal[i][botC12E2Index[index]].z); + if (DistVec2 <= dist + 0.5 && DistVec2 >= dist - 0.5) { + phiCount[dist].botPhiC12E2Count += 1; + } + } + + for (unsigned int index = 0; index < topC12E2MIndex.size(); index++) { + DistVec3 = trueDist(&NPX, &NPY, &NPZ, + &inputTotal[i][topC12E2MIndex[index]].x, + &inputTotal[i][topC12E2MIndex[index]].y, + &inputTotal[i][topC12E2MIndex[index]].z); + if (DistVec3 <= dist + 0.5 && DistVec3 >= dist - 0.5) { + phiCount[dist].topPhiC12E2MCount += 1; + } + } + + for (unsigned int index = 0; index < botC12E2MIndex.size(); index++) { + DistVec4 = trueDist(&NPX, &NPY, &NPZ, + &inputTotal[i][botC12E2MIndex[index]].x, + &inputTotal[i][botC12E2MIndex[index]].y, + &inputTotal[i][botC12E2MIndex[index]].z); + if (DistVec4 <= dist + 0.5 && DistVec4 >= dist - 0.5) { + phiCount[dist].botPhiC12E2MCount += 1; + } + } + } + phiTotal.push_back(phiCount); + phiCount.clear(); + } + } + + void PhiPrint() { + phipm A; + for (unsigned phiIndex = 0; phiIndex <= 100; ++phiIndex) { + A.phim = 0.0; + A.phip = 0.0; + A.phimVec.clear(); + A.phipVec.clear(); + NewNew.push_back(A); + } + + double phi1, phi2; + double phip, phim; + + for (unsigned int index = 0; index <= phiTotal.size() - 1; index++) { + for (unsigned int index2 = 0; index2 <= phiTotal[0].size(); index2++) { + + phi1 = 0.0; + phi2 = 0.0; + phip = 0.0; + phim = 0.0; + + // std::cout << index << " " << index2 << " " << + // phiTotal[index][index2].topPhiC12E2Count << " " << + // phiTotal[index][index2].botPhiC12E2Count << " " << + // phiTotal[index][index2].topPhiC12E2MCount << " " << + // phiTotal[index][index2].botPhiC12E2MCount << std::endl; + if ((double)(phiTotal[index][index2].topPhiC12E2Count + + phiTotal[index][index2].topPhiC12E2MCount) == 0.0 || + (double)(phiTotal[index][index2].botPhiC12E2Count + + phiTotal[index][index2].botPhiC12E2MCount) == 0) { + // Do nothing + + } else { + phi1 = (double)(phiTotal[index][index2].topPhiC12E2Count - + phiTotal[index][index2].topPhiC12E2MCount) / + (double)(phiTotal[index][index2].topPhiC12E2Count + + phiTotal[index][index2].topPhiC12E2MCount); + phi2 = (double)(phiTotal[index][index2].botPhiC12E2Count - + phiTotal[index][index2].botPhiC12E2MCount) / + (double)(phiTotal[index][index2].botPhiC12E2Count + + phiTotal[index][index2].botPhiC12E2MCount); + } + + phip = phi2 + phi1 / 2.0; + phim = phi2 - phi1 / 2.0; + + NewNew[index2].phim += phim; + NewNew[index2].phip += phip; + NewNew[index2].phimVec.push_back(phim); + NewNew[index2].phipVec.push_back(phip); + + std::cout << index << " " << index2 << " " << phip << " " << phim + << std::endl; + } + } + } + + void ComputePhiStandardDev() { + double sum, mean, sq_sum, stdev; + for (std::vector::iterator it = NewNew.begin(); it != NewNew.end(); + it++) { + // std::cout << " " << it - NewNew.begin() << " " << + // (it->phim)/inputTotal.size() << " " << (it->phip)/inputTotal.size() << + // " " << std::endl; std::cout << (it->phimVec).size() << " " << + // (it->phipVec).size() << std::endl; + for (unsigned int i = 0; i != NewNew[it - NewNew.begin()].phimVec.size(); + ++i) { + + sum = std::accumulate(NewNew[it - NewNew.begin()].phimVec.begin(), + NewNew[it - NewNew.begin()].phimVec.end(), + 0.0); // Compute sum + + mean = sum / NewNew[it - NewNew.begin()].phimVec.size(); // Compute Mean + + sq_sum = std::inner_product(NewNew[it - NewNew.begin()].phimVec.begin(), + NewNew[it - NewNew.begin()].phimVec.end(), + NewNew[it - NewNew.begin()].phimVec.begin(), + 0.0); // Compute square sum + + stdev = std::sqrt(sq_sum / NewNew[it - NewNew.begin()].phimVec.size() - + mean * mean) / + (pow(numberofSS, 0.5)); + } + + std::cout << " " << it - NewNew.begin() << " " + << (it->phip) / inputTotal.size() << " " << stdev << std::endl; + } + } + + void ComputeOrderphobic() { // Computes the phi, or the mismatch between the + // bilayer leaflets around the NP + + double DIST, DIST2; // TODO + double DISTM, DISTM2; // TODO + + for (unsigned int i = 0; i < inputTotal.size(); ++i) { + + for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { + + for (unsigned int newindex = 0; newindex < C12E2IndexVector.size(); + ++newindex) { + + DIST = trueDist(&inputTotal[i][C12E2IndexVector[index] + 4].x, + &inputTotal[i][C12E2IndexVector[index] + 4].y, + &inputTotal[i][C12E2IndexVector[index] + 4].z, + &inputTotal[i][C12E2IndexVector[newindex] + 4].x, + &inputTotal[i][C12E2IndexVector[newindex] + 4].y, + &inputTotal[i][C12E2IndexVector[newindex] + 4].z); + + C12E2sample.index = C12E2IndexVector[index]; + C12E2sample.dist = DIST; + C12E2sample.selfXcoord = inputTotal[i][C12E2IndexVector[index] + 4].x; + C12E2sample.selfYcoord = inputTotal[i][C12E2IndexVector[index] + 4].y; + C12E2sample.selfZcoord = inputTotal[i][C12E2IndexVector[index] + 4].z; + C12E2sample.Xcoord = inputTotal[i][C12E2IndexVector[newindex] + 4].x; + C12E2sample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].y; + C12E2sample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].z; + C12E2orderphobic.push_back(C12E2sample); + } + + for (unsigned int newindex = 0; newindex < C12E2IndexVector.size(); + ++newindex) { + + DISTM = trueDist(&inputTotal[i][C12E2IndexVector[index] + 4].x, + &inputTotal[i][C12E2IndexVector[index] + 4].y, + &inputTotal[i][C12E2IndexVector[index] + 4].z, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].x, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].y, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].z); + + C12E2sample.index = C12E2IndexVector[index]; + C12E2sample.dist = DISTM; + C12E2sample.selfXcoord = inputTotal[i][C12E2IndexVector[index] + 4].x; + C12E2sample.selfYcoord = inputTotal[i][C12E2IndexVector[index] + 4].y; + C12E2sample.selfZcoord = inputTotal[i][C12E2IndexVector[index] + 4].z; + C12E2sample.Xcoord = inputTotal[i][C12E2MIndexVector[newindex] + 4].x; + C12E2sample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex] + 4].y; + C12E2sample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex] + 4].z; + C12E2orderphobic.push_back(C12E2sample); + } + + sort(C12E2orderphobic.begin(), C12E2orderphobic.end(), + compareByIndexOPh); + C12E2orderphobic.erase(C12E2orderphobic.begin()); + C12E2orderphobic.erase(C12E2orderphobic.begin() + 6, + C12E2orderphobic.end()); + // C12E2orderphobic.resize(6); + C12E2orderphobicVec.push_back(C12E2orderphobic); + + for (unsigned int newindex = 0; newindex < C12E2MIndexVector.size(); + ++newindex) { + DISTM2 = trueDist(&inputTotal[i][C12E2MIndexVector[index] + 4].x, + &inputTotal[i][C12E2MIndexVector[index] + 4].y, + &inputTotal[i][C12E2MIndexVector[index] + 4].z, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].x, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].y, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].z); + + C12E2Msample.index = C12E2MIndexVector[index]; + C12E2Msample.dist = DISTM2; + + C12E2Msample.selfXcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].x; + C12E2Msample.selfYcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].y; + C12E2Msample.selfZcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].z; + + C12E2Msample.Xcoord = + inputTotal[i][C12E2MIndexVector[newindex] + 4].x; + C12E2Msample.Ycoord = + inputTotal[i][C12E2MIndexVector[newindex] + 4].y; + C12E2Msample.Ycoord = + inputTotal[i][C12E2MIndexVector[newindex] + 4].z; + + C12E2Morderphobic.push_back(C12E2Msample); + } + + for (unsigned int newindex = 0; newindex < C12E2MIndexVector.size(); + ++newindex) { + DIST2 = trueDist(&inputTotal[i][C12E2MIndexVector[index] + 4].x, + &inputTotal[i][C12E2MIndexVector[index] + 4].y, + &inputTotal[i][C12E2MIndexVector[index] + 4].z, + &inputTotal[i][C12E2IndexVector[newindex] + 4].x, + &inputTotal[i][C12E2IndexVector[newindex] + 4].y, + &inputTotal[i][C12E2IndexVector[newindex] + 4].z); + + C12E2Msample.index = C12E2MIndexVector[index]; + + C12E2Msample.dist = DIST2; + + C12E2Msample.selfXcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].x; + C12E2Msample.selfYcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].y; + C12E2Msample.selfZcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].z; + + C12E2Msample.Xcoord = inputTotal[i][C12E2IndexVector[newindex] + 4].x; + C12E2Msample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].y; + C12E2Msample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].z; + + C12E2Morderphobic.push_back(C12E2Msample); + } + + // std::sort(c12E2Morderphobic.begin(), c12E2Morderphobic.end(), + // compareByIndexOPh); std::sort(C12E2Morderphobic.begin(), + // C12E2Morderphobic.end(), std::greater()); + // C12E2Morderphobic.erase(C12E2Morderphobic.begin()); + // C12E2Morderphobic.resize(6); + + sort(C12E2Morderphobic.begin(), C12E2Morderphobic.end(), + compareByIndexOPh); + C12E2Morderphobic.erase(C12E2Morderphobic.begin()); + C12E2Morderphobic.erase(C12E2Morderphobic.begin() + 6, + C12E2Morderphobic.end()); + + // C12E2Morderphobic.resize(6); + + C12E2MorderphobicVec.push_back(C12E2Morderphobic); + + C12E2orderphobic.clear(); + C12E2Morderphobic.clear(); + } + + // Final push_back + orderphobicC12E2.push_back(C12E2orderphobicVec); + orderphobicC12E2M.push_back(C12E2MorderphobicVec); + + C12E2MorderphobicVec.clear(); + C12E2orderphobicVec.clear(); + } + } + + void OrderphobicSort() { + // Computes the phi, or the mismatch between the + // bilayer leaflets around the NP + + OPVal Val; + std::vector tempVal; + double output; + + for (unsigned int i = 0; i < orderphobicC12E2.size() - 1; ++i) { + for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { + output = calcAngle(&orderphobicC12E2[i][index]); + + std::cout << output << " " << i << " " << index << " " + << orderphobicC12E2[i][index][0].selfXcoord << " " + << orderphobicC12E2[i][index][0].selfYcoord << " " + << orderphobicC12E2[i][index][0].selfZcoord << std::endl; + + Val.X = orderphobicC12E2[i][index][0].selfXcoord; + Val.Y = orderphobicC12E2[i][index][0].selfYcoord; + Val.Z = orderphobicC12E2[i][index][0].selfZcoord; + Val.orderphobicVal = output; + tempVal.push_back(Val); + } + + for (unsigned int index = 0; index < C12E2MIndexVector.size(); ++index) { + output = calcAngle(&orderphobicC12E2M[i][index]); + std::cout << output << " " << i << " " << index << " " + << orderphobicC12E2M[i][index][0].selfXcoord << " " + << orderphobicC12E2M[i][index][0].selfYcoord << " " + << orderphobicC12E2M[i][index][0].selfZcoord << std::endl; + Val.X = orderphobicC12E2M[i][index][0].selfXcoord; + Val.Y = orderphobicC12E2M[i][index][0].selfYcoord; + Val.Z = orderphobicC12E2M[i][index][0].selfZcoord; + Val.orderphobicVal = output; + tempVal.push_back(Val); + } + + orderphobicVectorFinal.push_back(tempVal); + tempVal.clear(); + } + } + + void printop() { + + OPVal PlaceHolder; + double distP; + for (unsigned int dist = 0; dist <= 100; ++dist) { + PlaceHolder.X = 0.0; + PlaceHolder.Y = 0.0; + PlaceHolder.Z = 0.0; + ABC.push_back(PlaceHolder); + } + + for (unsigned int i = 0; i != orderphobicVectorFinal.size(); i++) { + for (unsigned int j = 0; j != orderphobicVectorFinal[i].size(); j++) { + + NPX = inputTotal[i][71312].x; // x coordinate of the NP + NPY = inputTotal[i][71312].y; // y coordinate of the NP + NPZ = inputTotal[i][71312].z; // z coordinate of the NP + + // std::cout << i << " " << j << " " << orderphobicVectorFinal[i][j].X + // << " " << orderphobicVectorFinal[i][j].Y << " " << + // orderphobicVectorFinal[i][j].Z << " " << + // orderphobicVectorFinal[i][j].orderphobicVal << " " << NPX << " " << + // NPY << " " << NPZ << std::endl; + + distP = trueDist(&NPX, &NPY, &NPZ, &orderphobicVectorFinal[i][j].X, + &orderphobicVectorFinal[i][j].Y, + &orderphobicVectorFinal[i][j].Z); + + for (unsigned int dist = 0; dist <= 100; ++dist) { + if (distP >= dist - 0.5 && distP <= dist + 0.5) { + ABC[dist].VV.push_back(orderphobicVectorFinal[i][j].orderphobicVal); + } + } + } + } + } + + void LargePrint() { + double sum, mean; + + for (unsigned int i = 0; i != ABC.size(); i++) { + sum = std::accumulate(ABC[i].VV.begin(), ABC[i].VV.end(), + 0.0); // Compute sum + mean = sum / ABC[i].VV.size(); // Compute Mean + std::cout << i << " " << mean << " " << std::endl; + } + } + + void allocateCOM() { + CenterOfMass(&C12E2IndexVector, &C12E2MIndexVector, &inputTotal, &C12E2COM, + &C12E2MCOM, &C12E2TotalCOMArray, &C12E2MTotalCOMArray); + } + +private: + std::vector inputVectorStruct; // push back all structs + std::vector>> + orderphobicC12E2; // This needs to be a simplified + std::vector>> + orderphobicC12E2M; // This needs to be simplfied also.. + // Vectors to store trajectory values + std::vector inputVector; // push back all structs + std::vector> + inputTotal; // push back vector of structs + std::vector + C12E2IndexVector; // push back C12E2 bead 7 indices (headgroups) + std::vector + C12E2MIndexVector; // push back all C12E2M bead 13 indices (headgroups) + // Vectors to store the Centers of Mass + std::vector C12E2COM; + std::vector C12E2MCOM; + std::vector phiCount; + std::vector> phiTotal; + std::vector> C12E2TotalCOMArray; + std::vector> C12E2MTotalCOMArray; + // top head groups C12E2 + std::vector topC12E2Index; + // top head groups C12E2M + std::vector botC12E2Index; + // bot head groups C12E2M + std::vector topC12E2MIndex; + // bot head groups C12E2M + std::vector botC12E2MIndex; + // TODO + std::vector NewNew; // TODO - need to rename + + FILE *ipf; /* input file */ + int atomtype; + int index, l, n; + int nlines = numberofatoms + 9; + double x, y, z; /*coordinates for the atoms in the box*/ + double box1; + double box2; + double NPX, NPY, NPZ; + double boxlength[boxdim]; + char line[100]; + inputCoord inputline; + double DistVec1, DistVec2, DistVec3, DistVec4; + + OPHstruct C12E2sample; + OPHstruct C12E2Msample; + + std::vector> C12E2orderphobicVec; // TODO + std::vector> C12E2MorderphobicVec; // TODO + std::vector C12E2orderphobic; + std::vector C12E2Morderphobic; + std::vector> orderphobicVectorFinal; // TODO + std::vector ABC; +}; + +class testClass {}; + +// Intitiate test class + +compute C12E2PhiOrderphobic; + +int main(int argc, char *argv[]) { + + C12E2PhiOrderphobic.storeFile(); + C12E2PhiOrderphobic.sortVectors(); + C12E2PhiOrderphobic.check(); + C12E2PhiOrderphobic.headGroupVectorFormation(); + C12E2PhiOrderphobic.allocateCOM(); + C12E2PhiOrderphobic.ComputePhi(); + C12E2PhiOrderphobic.PhiPrint(); + C12E2PhiOrderphobic.ComputePhiStandardDev(); + C12E2PhiOrderphobic.ComputeOrderphobic(); + C12E2PhiOrderphobic.OrderphobicSort(); + C12E2PhiOrderphobic.printop(); + C12E2PhiOrderphobic.LargePrint(); + + return 0; +} diff --git a/C12E2_work/.clang-format b/C12E2_work/.clang-format new file mode 100644 index 00000000..d1f38edc --- /dev/null +++ b/C12E2_work/.clang-format @@ -0,0 +1,216 @@ +--- +Language: Cpp +# BasedOnStyle: LLVM +AccessModifierOffset: -2 +AlignAfterOpenBracket: Align +AlignArrayOfStructures: None +AlignConsecutiveAssignments: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + PadOperators: true +AlignConsecutiveBitFields: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + PadOperators: false +AlignConsecutiveDeclarations: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + PadOperators: false +AlignConsecutiveMacros: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + PadOperators: false +AlignEscapedNewlines: Right +AlignOperands: Align +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortEnumsOnASingleLine: true +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: Never +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: MultiLine +AttributeMacros: + - __capability +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Never + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: false + BeforeElse: false + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeConceptDeclarations: Always +BreakBeforeBraces: Attach +BreakBeforeInheritanceComma: false +BreakInheritanceList: BeforeColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 80 +CommentPragmas: '^ IWYU pragma:' +QualifierAlignment: Leave +CompactNamespaces: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DeriveLineEnding: true +DerivePointerAlignment: false +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +PackConstructorInitializers: BinPack +BasedOnStyle: '' +ConstructorInitializerAllOnOneLineOrOnePerLine: false +AllowAllConstructorInitializersOnNextLine: true +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Preserve +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + SortPriority: 0 + CaseSensitive: false + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + SortPriority: 0 + CaseSensitive: false + - Regex: '.*' + Priority: 1 + SortPriority: 0 + CaseSensitive: false +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseLabels: false +IndentCaseBlocks: false +IndentGotoLabels: true +IndentPPDirectives: None +IndentExternBlock: AfterExternBlock +IndentRequiresClause: true +IndentWidth: 2 +IndentWrappedFunctionNames: false +InsertBraces: false +InsertTrailingCommas: None +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +LambdaBodyIndentation: Signature +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 2 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakOpenParenthesis: 0 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PenaltyIndentedWhitespace: 0 +PointerAlignment: Right +PPIndentWidth: -1 +ReferenceAlignment: Pointer +ReflowComments: true +RemoveBracesLLVM: false +RequiresClausePosition: OwnLine +SeparateDefinitionBlocks: Leave +ShortNamespaceLines: 1 +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceBeforeParensOptions: + AfterControlStatements: true + AfterForeachMacros: true + AfterFunctionDefinitionName: false + AfterFunctionDeclarationName: false + AfterIfMacros: true + AfterOverloadedOperator: false + AfterRequiresInClause: false + AfterRequiresInExpression: false + BeforeNonEmptyParentheses: false +SpaceAroundPointerQualifiers: Default +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyBlock: false +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: Never +SpacesInConditionalStatement: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParentheses: false +SpacesInSquareBrackets: false +SpaceBeforeSquareBrackets: false +BitFieldColonSpacing: Both +Standard: Latest +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 8 +UseCRLF: false +UseTab: Never +WhitespaceSensitiveMacros: + - STRINGIZE + - PP_STRINGIZE + - BOOST_PP_STRINGIZE + - NS_SWIFT_NAME + - CF_SWIFT_NAME +... + diff --git a/C12E2_work/C12E2FullAnalysis.cxx b/C12E2_work/C12E2FullAnalysis.cxx index 5cb8d430..f3915efb 100755 --- a/C12E2_work/C12E2FullAnalysis.cxx +++ b/C12E2_work/C12E2FullAnalysis.cxx @@ -1,60 +1,62 @@ /* - ------------------------------------------------------------------------ - | Analysis of C12E2 Mixed Bilayers | ------------------------------------------------------------------------ - + | Analysis of C12E2 Mixed Bilayers | + ------------------------------------------------------------------------ + Author = "Sang Young Noh" Version = "0.0.1" - - Updated = 02/03/2019 + + Updated = 11/09/2019 The papers referenced for this work is: - 1. "F. C. MacKintosh and S. A. Safran, Phase Separation and curvature of bilayer membranes, Physical Review E, 47, 2, 1993" + 1. "F. C. MacKintosh and S. A. Safran, Phase Separation and curvature of + bilayer membranes, Physical Review E, 47, 2, 1993" - 2. "R. Lipowsky, Domain-induced budding of fluid membranes, Biophysical Society, 64, 1993, 1133 - 1138" + 2. "R. Lipowsky, Domain-induced budding of fluid membranes, Biophysical + Society, 64, 1993, 1133 - 1138" - 3. "J. Wolff and S. Komura and D. Andelman, Budding of domains in mixed bilayer domains, 91, Physical Review E, 91, 012708, 2015" - - 4. S. Katira and K. K. Madadapu and S. Vaikuntanathan and B. Smit and D. Chandler, Pre-transition effects mediate forces of assembly between transmembrane proteins, elife, 2016, 5, e13150 + 3. "J. Wolff and S. Komura and D. Andelman, Budding of domains in mixed bilayer + domains, 91, Physical Review E, 91, 012708, 2015" - 5. Physical Considerations of the Organization of Inclusions in Lipid Bilayer Systems, Spring 2015, Shachi Katira + 4. S. Katira and K. K. Madadapu and S. Vaikuntanathan and B. Smit and D. + Chandler, Pre-transition effects mediate forces of assembly between + transmembrane proteins, elife, 2016, 5, e13150 -*/ + 5. Physical Considerations of the Organization of Inclusions in Lipid Bilayer + Systems, Spring 2015, Shachi Katira -#include -#include +*/ -#include #include -#include // std::multiplies -#include // std::adjacent_differenc#include -#include -#include +#include +#include #include -#include #include -#include -#include -#include - #include // Library for complex numbers +#include +#include // std::multiplies +#include +#include +#include +#include // std::adjacent_differenc#include +#include +#include +#include #define _USE_MATH_DEFINES -//#include -//#include -//using namespace Eigen; - -const int numberOfPolymers = 998; // The number of polymers of each type - C12E2 or mimic +// Constants for use later with analysis +const int NumberOfPolymers = + 998; // The number of polymers of each type - C12E2 or mimic const int numberofatoms = 71313; // Total number of beads in the simulation const int indexCG = 7; int numberofSS = 100; /*The number of screenshots in the dump file*/ const int boxdim = 3; -typedef struct { +typedef struct { int index[7]; } C12E2_skeleton; @@ -73,7 +75,7 @@ typedef struct { // Used to input the center of masses for each lipid double z; } COMstruct; -typedef struct { // Used to input the center of masses for each lipid +typedef struct { double dist; int topPhiC12E2Count; int topPhiC12E2MCount; @@ -81,21 +83,21 @@ typedef struct { // Used to input the center of masses for each lipid int botPhiC12E2MCount; } phiStruct; -typedef struct { // Used to input the center of masses for each lipid +typedef struct { double phip; double phim; std::vector phipVec; std::vector phimVec; } phipm; -typedef struct { // Used to input the center of masses for each lipid +typedef struct { double X, Y, Z; double orderphobicVal; - std::vector VV; + std::vector surrounding_residue_indices; } OPVal; - -typedef struct { // Used to identify the group and distance to compute the orderphobic effect +typedef struct { // Used to identify the group and distance to compute the + // orderphobic effect int index; double dist; double selfXcoord; @@ -104,391 +106,491 @@ typedef struct { // Used to identify the group and distance to compute the order double Xcoord; double Ycoord; double Zcoord; - + } OPHstruct; +double calcAngle(std::vector *OPHInput) { -/* -double CalculateOrderphobicEffect() { - std::complex - - return double; -} -*/ -double calcAngle(std::vector* OPHInput) { - double refVector[2]; - double dot; + double dot; double newx; double newy; double u, v; + std::complex angle; - double phi; + + double phi; + refVector[0] = 10.0; + refVector[1] = 10.0; - std::complex mycomplex(0.0, 0.0); - std::complex sixth((double)1/6, 0.0); + + std::complex mycomplex(0.0, 0.0); + std::complex sixth((double)1 / 6, 0.0); + double absval; + for (unsigned int i = 0; i < OPHInput->size(); ++i) { + newx = OPHInput->at(i).Xcoord; + newy = OPHInput->at(i).Ycoord; - u = pow((pow(refVector[0],2) + pow(refVector[1],2)), 0.5); - v = pow((pow(newx,2) + pow(newy,2)), 0.5); + + u = pow((pow(refVector[0], 2) + pow(refVector[1], 2)), 0.5); + + v = pow((pow(newx, 2) + pow(newy, 2)), 0.5); + dot = (newx * refVector[0]) + (newy * refVector[1]); - std::complex placeholder(0, ( 6 * 180 * M_PI * acos((dot / (u * v))))); // TODO - angle = exp(placeholder); + + std::complex placeholder( + 0, (6 * 180 * M_PI * acos((dot / (u * v))))); // TODO + + angle = exp(placeholder); mycomplex += angle; - // std::cout << i << " " << absval << std::endl; + + // std::cout << i << " " << absval << std::endl; } mycomplex = sixth * mycomplex; - absval = pow(abs(mycomplex),2); - //std::cout << absval << std::endl; + absval = pow(abs(mycomplex), 2); + // std::cout << absval << std::endl; return absval; } -double trueDist(double* COM1x, double* COM1y, double* COM1z, double* COM2x, double* COM2y, double* COM2z) { - double dist = pow((pow(*COM1x-*COM2x,2) + pow(*COM1y-*COM2y,2) + pow(*COM1z-*COM2z,2)),0.5); +double trueDist(double *COM1x, double *COM1y, double *COM1z, double *COM2x, + double *COM2y, double *COM2z) { + double dist = pow((pow(*COM1x - *COM2x, 2) + pow(*COM1y - *COM2y, 2) + + pow(*COM1z - *COM2z, 2)), + 0.5); return dist; } -// Function used for sorting vector of structs +// Function used for sorting vector of structs bool compareByIndex(const inputCoord &a, const inputCoord &b) { - return a.a < b.a; + return a.a < b.a; } -bool compareByIndexOPh(const OPHstruct &a, const OPHstruct &b) { - return a.dist < b.dist; +bool compareByIndexOPh(const OPHstruct &a, const OPHstruct &b) { + return a.dist < b.dist; } +void CenterOfMass(std::vector *vec1, std::vector *vec2, + std::vector> *inputVec, + std::vector *COM1, std::vector *COM2, + std::vector> *C12E2final, + std::vector> *C12E2Mfinal) { + // CenterOfMass(&C12E2IndexVector, &C12E2MIndexVector, &inputTotal, + // &C12E2COM, &C12E2MCOM, &C12E2TotalCOMArray, &C12E2MTotalCOMArray); -void CenterOfMass(std::vector* vec1, std::vector* vec2, std::vector >* inputVec, std::vector* COM1, std::vector* COM2, std::vector >* C12E2final, std::vector >* C12E2Mfinal) { - // CenterOfMass(&C12E2IndexVector, &C12E2MIndexVector, &inputTotal, &C12E2COM, &C12E2MCOM, &C12E2TotalCOMArray, &C12E2MTotalCOMArray); - - double C12E2comX, C12E2McomX; // X coordinate COM - double C12E2comY, C12E2McomY; // Y coordinate COM + double C12E2comX, C12E2McomX; // X coordinate COM + double C12E2comY, C12E2McomY; // Y coordinate COM double C12E2comZ, C12E2McomZ; // Z coordinate cCOM - - COMstruct C12E2input; // struct to store COM coordinates for C12E2 - COMstruct C12E2Minput; // struct to store cOM coordinates for C12E2M - - for (unsigned int i = 0; i <= inputVec->size()-1; ++i) { - - for (unsigned int index = 0; index <= vec1->size()-1; ++index) { - - C12E2comX = (inputVec->at(i).at(vec1->at(index)).x + inputVec->at(i).at((vec1->at(index))+1).x + inputVec->at(i).at((vec1->at(index))+2).x + inputVec->at(i).at((vec1->at(index))+3).x + inputVec->at(i).at((vec1->at(index))+4).x + inputVec->at(i).at((vec1->at(index))+5).x + inputVec->at(i).at((vec1->at(index))+6).x)/7.0; - C12E2comY = (inputVec->at(i).at(vec1->at(index)).y + inputVec->at(i).at((vec1->at(index))+1).x + inputVec->at(i).at((vec1->at(index))+2).y + inputVec->at(i).at((vec1->at(index))+3).y + inputVec->at(i).at((vec1->at(index))+4).y + inputVec->at(i).at((vec1->at(index))+5).y + inputVec->at(i).at((vec1->at(index))+6).y)/7.0; + COMstruct C12E2input; // struct to store COM coordinates for C12E2 + COMstruct C12E2Minput; // struct to store cOM coordinates for C12E2M - C12E2comZ = (inputVec->at(i).at(vec1->at(index)).z + inputVec->at(i).at((vec1->at(index))+1).z + inputVec->at(i).at((vec1->at(index))+2).z + inputVec->at(i).at((vec1->at(index))+3).z + inputVec->at(i).at((vec1->at(index))+4).z + inputVec->at(i).at((vec1->at(index))+5).z + inputVec->at(i).at((vec1->at(index))+6).z)/7.0; + for (unsigned int i = 0; i <= inputVec->size() - 1; ++i) { + + for (unsigned int index = 0; index <= vec1->size() - 1; ++index) { + + C12E2comX = (inputVec->at(i).at(vec1->at(index)).x + + inputVec->at(i).at((vec1->at(index)) + 1).x + + inputVec->at(i).at((vec1->at(index)) + 2).x + + inputVec->at(i).at((vec1->at(index)) + 3).x + + inputVec->at(i).at((vec1->at(index)) + 4).x + + inputVec->at(i).at((vec1->at(index)) + 5).x + + inputVec->at(i).at((vec1->at(index)) + 6).x) / + 7.0; + + C12E2comY = (inputVec->at(i).at(vec1->at(index)).y + + inputVec->at(i).at((vec1->at(index)) + 1).x + + inputVec->at(i).at((vec1->at(index)) + 2).y + + inputVec->at(i).at((vec1->at(index)) + 3).y + + inputVec->at(i).at((vec1->at(index)) + 4).y + + inputVec->at(i).at((vec1->at(index)) + 5).y + + inputVec->at(i).at((vec1->at(index)) + 6).y) / + 7.0; + + C12E2comZ = (inputVec->at(i).at(vec1->at(index)).z + + inputVec->at(i).at((vec1->at(index)) + 1).z + + inputVec->at(i).at((vec1->at(index)) + 2).z + + inputVec->at(i).at((vec1->at(index)) + 3).z + + inputVec->at(i).at((vec1->at(index)) + 4).z + + inputVec->at(i).at((vec1->at(index)) + 5).z + + inputVec->at(i).at((vec1->at(index)) + 6).z) / + 7.0; C12E2input.index = vec1->at(index); C12E2input.x = C12E2comX; C12E2input.y = C12E2comY; C12E2input.z = C12E2comZ; - COM1->push_back(C12E2input); + COM1->push_back(C12E2input); } - - for (unsigned int index = 0; index <= vec2->size()-1; ++index) { - - C12E2McomX = (inputVec->at(i).at(vec2->at(index)).x + inputVec->at(i).at((vec2->at(index))+1).x + inputVec->at(i).at((vec2->at(index))+2).x + inputVec->at(i).at((vec2->at(index))+3).x + inputVec->at(i).at((vec2->at(index))+4).x + inputVec->at(i).at((vec2->at(index))+5).x + inputVec->at(i).at((vec2->at(index))+6).x)/7.0; - - C12E2McomY = (inputVec->at(i).at(vec2->at(index)).y + inputVec->at(i).at((vec2->at(index))+1).x + inputVec->at(i).at((vec2->at(index))+2).y + inputVec->at(i).at((vec2->at(index))+3).y + inputVec->at(i).at((vec2->at(index))+4).y + inputVec->at(i).at((vec2->at(index))+5).y + inputVec->at(i).at((vec2->at(index))+6).y)/7.0; - C12E2McomZ = (inputVec->at(i).at(vec2->at(index)).z + inputVec->at(i).at((vec2->at(index))+1).z + inputVec->at(i).at((vec2->at(index))+2).z + inputVec->at(i).at((vec2->at(index))+3).z + inputVec->at(i).at((vec2->at(index))+4).z + inputVec->at(i).at((vec2->at(index))+5).z + inputVec->at(i).at((vec2->at(index))+6).z)/7.0; + for (unsigned int index = 0; index <= vec2->size() - 1; ++index) { + + C12E2McomX = (inputVec->at(i).at(vec2->at(index)).x + + inputVec->at(i).at((vec2->at(index)) + 1).x + + inputVec->at(i).at((vec2->at(index)) + 2).x + + inputVec->at(i).at((vec2->at(index)) + 3).x + + inputVec->at(i).at((vec2->at(index)) + 4).x + + inputVec->at(i).at((vec2->at(index)) + 5).x + + inputVec->at(i).at((vec2->at(index)) + 6).x) / + 7.0; + + C12E2McomY = (inputVec->at(i).at(vec2->at(index)).y + + inputVec->at(i).at((vec2->at(index)) + 1).x + + inputVec->at(i).at((vec2->at(index)) + 2).y + + inputVec->at(i).at((vec2->at(index)) + 3).y + + inputVec->at(i).at((vec2->at(index)) + 4).y + + inputVec->at(i).at((vec2->at(index)) + 5).y + + inputVec->at(i).at((vec2->at(index)) + 6).y) / + 7.0; + + C12E2McomZ = (inputVec->at(i).at(vec2->at(index)).z + + inputVec->at(i).at((vec2->at(index)) + 1).z + + inputVec->at(i).at((vec2->at(index)) + 2).z + + inputVec->at(i).at((vec2->at(index)) + 3).z + + inputVec->at(i).at((vec2->at(index)) + 4).z + + inputVec->at(i).at((vec2->at(index)) + 5).z + + inputVec->at(i).at((vec2->at(index)) + 6).z) / + 7.0; C12E2Minput.index = vec2->at(index); C12E2Minput.x = C12E2comX; C12E2Minput.y = C12E2comY; C12E2Minput.z = C12E2comZ; - //TODO - COM1->push_back(C12E2Minput); + // TODO + COM1->push_back(C12E2Minput); } - + C12E2final->push_back(*COM1); // C12E2final->push_back(*COM2); - //TODO + // TODO C12E2Mfinal->push_back(*COM2); COM1->clear(); COM2->clear(); } } -bool sortbysec_int(const std::pair &a, const std::pair &b) { - return (a.second < b.second); +bool sortbysec_int(const std::pair &a, const std::pair &b) { + return (a.second < b.second); } -bool sortbysec_double(const std::pair &a, const std::pair &b) { - return (a.second < b.second); -} +bool sortbysec_double(const std::pair &a, + const std::pair &b) { + return (a.second < b.second); +} -//pointCluster* regionQuery (int , double*, double*, double*, double*, double*, double*); -//#define minClust 5 // the minimum sizes cluster we want to count as a cluster class compute { public: - compute() {}; // Default constructor - compute(std::string file) {}; // constructor reading file - void storeFile() { + compute(){}; // Default constructor + compute(std::string file){}; // constructor reading file + + void storeFile() { + boost::progress_display show_progress(numberofSS); - // open file for reading - ipf = fopen("dump.mixed2", "r"); // Needs correction - // check for error opening file */ - if (ipf == NULL) { + // open file for reading + ipf = fopen("dump.mixed2", "r"); // Needs correction + // check for error opening file */ + if (ipf == NULL) { std::cout << "Error opening file\n"; exit(1); } - // loop over the values - for (int SSno = 0; SSno < numberofSS; SSno++) { + // loop over the values + for (int SSno = 0; SSno < numberofSS; SSno++) { l = 0; n = 0; index = 0; // Make sure to clear all vectors before going on to the next snapshot inputVector.clear(); - for (int k = 0; k < nlines; k++) { - // get a line from the file - // fgets() returns NULL when it reaches an error or end of file - - fgets(line,sizeof(line),ipf); - //while (fgets(line,sizeof(line),ipf) != NULL) { - - if (l < 5) { - /* We are doing nothing */ - // printf("%s l= %d", line,l); - l++; - } - else if ((l > 4 && l < 8)) { - /*We are scanning the bit with just the box parameters*/ - sscanf(line, "%lf %lf", &box1, &box2); - // printf("%lf %lf \n",box1,box2 ); - boxlength[l-5] = box2-box1; - l++; - } - else if (l == 8) { - // printf(" **** l= %d \n",l); - /* We are doing nothing */ - l++; - } - else { - // printf(" ***** l = %d \n",l ); - /* convert the text to numbers */ - - sscanf(line,"%d %d %lf %lf %lf",&index,&atomtype,&x,&y,&z); - //std::cout << index << " " << atomtype << " " << x << " " << y << " " << z; - //a.push_back(std::make_pair(index, index)); // Push back indices - //b.push_back(std::make_pair(index, atomtype)); // Push back atomtypes - //xco.push_back(std::make_pair(index, x*boxlength[0])); // Push back boxlengths - x coordinates - //yco.push_back(std::make_pair(index, y*boxlength[1])); // Push back boxlengths - y coordinates - //zco.push_back(std::make_pair(index, z*boxlength[2])); // Push back boxlengths - z coordinates - inputline.a = index; - inputline.b = atomtype; - inputline.x = x*boxlength[0]; - inputline.y = y*boxlength[1]; - inputline.z = z*boxlength[2]; - inputVector.push_back(inputline); - n++; - l++; - //printf("%d %d %lf %lf %lf\n",index,atomtype,x,y,z); - } - + for (int k = 0; k < nlines; k++) { + // get a line from the file + // fgets() returns NULL when it reaches an error or end of file + + fgets(line, sizeof(line), ipf); + // while (fgets(line,sizeof(line),ipf) != NULL) { + + if (l < 5) { + /* We are doing nothing */ + // printf("%s l= %d", line,l); + l++; + } else if ((l > 4 && l < 8)) { + /*We are scanning the bit with just the box parameters*/ + sscanf(line, "%lf %lf", &box1, &box2); + // printf("%lf %lf \n",box1,box2 ); + boxlength[l - 5] = box2 - box1; + l++; + } else if (l == 8) { + // printf(" **** l= %d \n",l); + /* We are doing nothing */ + l++; + } else { + // printf(" ***** l = %d \n",l ); + /* convert the text to numbers */ + + sscanf(line, "%d %d %lf %lf %lf", &index, &atomtype, &x, &y, &z); + // std::cout << index << " " << atomtype << " " << x << " " << y << " + // " << z; a.push_back(std::make_pair(index, index)); // Push back + // indices b.push_back(std::make_pair(index, atomtype)); // Push back + // atomtypes xco.push_back(std::make_pair(index, x*boxlength[0])); // + // Push back boxlengths - x coordinates + // yco.push_back(std::make_pair(index, y*boxlength[1])); // Push back + // boxlengths - y coordinates zco.push_back(std::make_pair(index, + // z*boxlength[2])); // Push back boxlengths - z coordinates + inputline.a = index; + inputline.b = atomtype; + inputline.x = x * boxlength[0]; + inputline.y = y * boxlength[1]; + inputline.z = z * boxlength[2]; + inputVector.push_back(inputline); + n++; + l++; + // printf("%d %d %lf %lf %lf\n",index,atomtype,x,y,z); + } } inputTotal.push_back(inputVector); ++show_progress; } } - void sortVectors () { // Sort vector values based on indices of the dump files + void sortVectors() { // Sort vector values based on indices of the dump files for (unsigned int i = 0; i < inputTotal.size(); ++i) { sort(inputTotal[i].begin(), inputTotal[i].end(), compareByIndex); } } - - void printVectorElements() { - for (unsigned int i = 0; i < inputTotal.size(); ++i) { - for (unsigned int j = 0; j < inputTotal[1].size(); ++j) { - //std::cout << i << " " << j << " " << " " << inputTotal[i][j].a << " " << inputTotal[i][j].b << " " << inputTotal[i][j].x << " " << inputTotal[i][j].y << " " << inputTotal[i][j].z << std::endl; - } - } + void printVectorElements() { + for (unsigned int i = 0; i < inputTotal.size(); ++i) { + for (unsigned int j = 0; j < inputTotal[1].size(); ++j) { + // std::cout << i << " " << j << " " << " " << inputTotal[i][j].a << " " + // << inputTotal[i][j].b << " " << inputTotal[i][j].x << " " << + // inputTotal[i][j].y << " " << inputTotal[i][j].z << std::endl; + } + } } - - /* - void CompositionProfile() { - //Calculating the composition profile - for (int i = 0; i <= sizeof(C12E2_struct)/sizeof(C12E2_struct[1])-1; i++) { // TODO - //double truedist(double COM1x, double COM1y, double COM1z, double COM2x, double COM2y, double COM2z) { - // See if the distance between the nanopparticle is smaller than 15 + /* + void CompositionProfile() { + //Calculating the composition profile + for (int i = 0; i <= sizeof(C12E2_struct)/sizeof(C12E2_struct[1])-1; i++) { + // TODO + + //double truedist(double COM1x, double COM1y, double COM1z, double COM2x, + double COM2y, double COM2z) { + // See if the distance between the nanopparticle is smaller than 15 //std::cout << zco[71312] << std::endl; - if (trueDist(xco[71312], yco[71312], zco[71312], headGroupC12_E2xCOM[i], headGroupC12_E2yCOM[i], headGroupC12_E2zCOM[i]) <= 25.000) { - // Check if the c12e2 molecule is on the top layer or the bottom layer - if (headGroupC12_E2zCOM[i] > zco[71312]) { // If on the top layer - tophead++; - } - else if (headGroupC12_E2zCOM[i] < zco[71312]) { // If on the bottom layer - downhead++; - } + if (trueDist(xco[71312], yco[71312], zco[71312], headGroupC12_E2xCOM[i], + headGroupC12_E2yCOM[i], headGroupC12_E2zCOM[i]) <= 25.000) { + // Check if the c12e2 molecule is on the top layer or the bottom layer + if (headGroupC12_E2zCOM[i] > zco[71312]) { // If on the top layer + tophead++; + } + else if (headGroupC12_E2zCOM[i] < zco[71312]) { // If on the bottom + layer downhead++; + } } - - // This time, working with mimics - if(trueDist(xco[71312], yco[71312], zco[71312], headGroupMIMCxCOM[i], headGroupMIMICyCOM[i], headGroupMIMICzCOM[i] ) <= 25.000) { - if (headGroupMIMICzCOM[i] > zco[71312]) { - mimictophead++; - } - else if (headGroupMIMICzCOM[i] < zco[71312]) { - mimicdownhead++; - } - } + + // This time, working with mimics + if(trueDist(xco[71312], yco[71312], zco[71312], headGroupMIMCxCOM[i], + headGroupMIMICyCOM[i], headGroupMIMICzCOM[i] ) <= 25.000) { if + (headGroupMIMICzCOM[i] > zco[71312]) { mimictophead++; + } + else if (headGroupMIMICzCOM[i] < zco[71312]) { + mimicdownhead++; + } + } } - - - //std::cout << SSno << " " << tophead/(tophead + mimictophead) << " " << mimictophead/(tophead + mimictophead) << " " << downhead/(downhead + mimicdownhead) << " " << mimicdownhead/(downhead + mimicdownhead) << " " << (tophead-mimictophead)/(tophead + mimictophead) << " " << (downhead-mimicdownhead)/(downhead + mimicdownhead) << " " << ((tophead-mimictophead)/(tophead + mimictophead) + (downhead-mimicdownhead)/(downhead + mimicdownhead))/2- << " " << ((tophead-mimictophead)/(tophead + mimictophead) - (downhead-mimicdownhead)/(downhead + mimicdownhead))/2 << std::endl; + + + //std::cout << SSno << " " << tophead/(tophead + mimictophead) << " " << + mimictophead/(tophead + mimictophead) << " " << downhead/(downhead + + mimicdownhead) << " " << mimicdownhead/(downhead + mimicdownhead) << " " << + (tophead-mimictophead)/(tophead + mimictophead) << " " << + (downhead-mimicdownhead)/(downhead + mimicdownhead) << " " << + ((tophead-mimictophead)/(tophead + mimictophead) + + (downhead-mimicdownhead)/(downhead + mimicdownhead))/2- << " " << + ((tophead-mimictophead)/(tophead + mimictophead) - + (downhead-mimicdownhead)/(downhead + mimicdownhead))/2 << std::endl; } */ - - void check() { - for (unsigned int i = 0; i < inputTotal.size(); ++i) { + + void check() { + for (unsigned int i = 0; i < inputTotal.size(); ++i) { for (unsigned int j = 0; j <= inputTotal[1].size(); j++) { - if (inputTotal[i][j].b == 7) { - std::cout << inputTotal[i][j].a << " " << inputTotal[i][j].b << std::endl; - std::cout << inputTotal[i][j+1].a << " " << inputTotal[i][j+1].b << std::endl; - std::cout << inputTotal[i][j+2].a << " " << inputTotal[i][j+2].b << std::endl; - std::cout << inputTotal[i][j+3].a << " " << inputTotal[i][j+3].b << std::endl; - std::cout << inputTotal[i][j+4].a << " " << inputTotal[i][j+4].b << std::endl; - std::cout << inputTotal[i][j+5].a << " " << inputTotal[i][j+5].b << std::endl; - std::cout << inputTotal[i][j+6].a << " " << inputTotal[i][j+6].b << std::endl; - } else if (inputTotal[i][j].b == 13) { - std::cout << inputTotal[i][j].a << " " << inputTotal[i][j].b << std::endl; - std::cout << inputTotal[i][j+1].a << " " << inputTotal[i][j+1].b << std::endl; - std::cout << inputTotal[i][j+2].a << " " << inputTotal[i][j+2].b << std::endl; - std::cout << inputTotal[i][j+3].a << " " << inputTotal[i][j+3].b << std::endl; - std::cout << inputTotal[i][j+4].a << " " << inputTotal[i][j+4].b << std::endl; - std::cout << inputTotal[i][j+5].a << " " << inputTotal[i][j+5].b << std::endl; - std::cout << inputTotal[i][j+6].a << " " << inputTotal[i][j+6].b << std::endl; - } + if (inputTotal[i][j].b == 7) { + std::cout << inputTotal[i][j].a << " " << inputTotal[i][j].b + << std::endl; + std::cout << inputTotal[i][j + 1].a << " " << inputTotal[i][j + 1].b + << std::endl; + std::cout << inputTotal[i][j + 2].a << " " << inputTotal[i][j + 2].b + << std::endl; + std::cout << inputTotal[i][j + 3].a << " " << inputTotal[i][j + 3].b + << std::endl; + std::cout << inputTotal[i][j + 4].a << " " << inputTotal[i][j + 4].b + << std::endl; + std::cout << inputTotal[i][j + 5].a << " " << inputTotal[i][j + 5].b + << std::endl; + std::cout << inputTotal[i][j + 6].a << " " << inputTotal[i][j + 6].b + << std::endl; + } else if (inputTotal[i][j].b == 13) { + std::cout << inputTotal[i][j].a << " " << inputTotal[i][j].b + << std::endl; + std::cout << inputTotal[i][j + 1].a << " " << inputTotal[i][j + 1].b + << std::endl; + std::cout << inputTotal[i][j + 2].a << " " << inputTotal[i][j + 2].b + << std::endl; + std::cout << inputTotal[i][j + 3].a << " " << inputTotal[i][j + 3].b + << std::endl; + std::cout << inputTotal[i][j + 4].a << " " << inputTotal[i][j + 4].b + << std::endl; + std::cout << inputTotal[i][j + 5].a << " " << inputTotal[i][j + 5].b + << std::endl; + std::cout << inputTotal[i][j + 6].a << " " << inputTotal[i][j + 6].b + << std::endl; + } } } } - void headGroupVectorFormation() { - for (unsigned int j = 0; j <= inputTotal[1].size(); j++) { + void headGroupVectorFormation() { + for (unsigned int j = 0; j <= inputTotal[1].size(); j++) { if (inputTotal[1][j].b == 7) { - C12E2IndexVector.push_back(inputTotal[1][j].a); // push back C12E2 bead 7 indices (headgroups) + C12E2IndexVector.push_back( + inputTotal[1][j].a); // push back C12E2 bead 7 indices (headgroups) } else if (inputTotal[1][j].b == 13) { - C12E2MIndexVector.push_back(inputTotal[1][j].a); // push back C12E2 bead 7 indices (headgroups) + C12E2MIndexVector.push_back( + inputTotal[1][j].a); // push back C12E2 bead 7 indices (headgroups) } } - + for (unsigned int index = 0; index <= C12E2IndexVector.size(); ++index) { - std::cout << C12E2IndexVector[index] << std::endl; + std::cout << C12E2IndexVector[index] << std::endl; } - for (unsigned int index = 0; index <= C12E2MIndexVector.size(); ++index) { - std::cout << C12E2MIndexVector[index] << std::endl; + for (unsigned int index = 0; index <= C12E2MIndexVector.size(); ++index) { + std::cout << C12E2MIndexVector[index] << std::endl; } } - - void ComputePhi() { // Computes the phi, or the mismatch between the bilayer leaflets around the NP - //std::cout << C12E2IndexVector.size() << " " << C12E2MIndexVector.size() << " " << inputTotal.size() << std::endl; + + void ComputePhi() { // Computes the phi, or the mismatch between the bilayer + // leaflets around the NP + // std::cout << C12E2IndexVector.size() << " " << C12E2MIndexVector.size() + // << " " << inputTotal.size() << std::endl; /* - Following Python's comment system """ """ + Following Python's comment system """ """ - We are explcitily ignoring flip-flops in the case of this study + We are explcitily ignoring flip-flops in the case of this study */ - - for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { - if (inputTotal[1][C12E2IndexVector[index]].z > 120.0) { - topC12E2Index.push_back(C12E2IndexVector[index]); - } else if (inputTotal[1][C12E2IndexVector[index]].z < 120.0) { - botC12E2Index.push_back(C12E2IndexVector[index]); + + for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { + if (inputTotal[1][C12E2IndexVector[index]].z > 120.0) { + + topC12E2Index.push_back(C12E2IndexVector[index]); + } else if (inputTotal[1][C12E2IndexVector[index]].z < 120.0) { + botC12E2Index.push_back(C12E2IndexVector[index]); } } - for (unsigned int index = 0; index < C12E2MIndexVector.size(); ++index) { - if (inputTotal[1][C12E2MIndexVector[index]].z > 120.0) { - topC12E2MIndex.push_back(C12E2MIndexVector[index]); - } else if (inputTotal[1][C12E2MIndexVector[index]].z < 120.0) { - botC12E2MIndex.push_back(C12E2MIndexVector[index]); - } + for (unsigned int index = 0; index < C12E2MIndexVector.size(); ++index) { + if (inputTotal[1][C12E2MIndexVector[index]].z > 120.0) { + topC12E2MIndex.push_back(C12E2MIndexVector[index]); + } else if (inputTotal[1][C12E2MIndexVector[index]].z < 120.0) { + botC12E2MIndex.push_back(C12E2MIndexVector[index]); } + } - std::cout << "check" << std::endl; + std::cout << "check" << std::endl; - for (std::vector::iterator it = topC12E2Index.begin(); it != topC12E2Index.end(); it++) { - //std::cout << *it << " " << std::endl; - // std::cout << topC12E2Index.size() << std::endl; - } + for (std::vector::iterator it = topC12E2Index.begin(); + it != topC12E2Index.end(); it++) { + // std::cout << *it << " " << std::endl; + // std::cout << topC12E2Index.size() << std::endl; + } - for (std::vector::iterator it = botC12E2Index.begin(); it != botC12E2Index.end(); it++) { - //std::cout << *it << " " << std::endl; - // std::cout << botC12E2Index.size() << std::endl; - } + for (std::vector::iterator it = botC12E2Index.begin(); + it != botC12E2Index.end(); it++) { + // std::cout << *it << " " << std::endl; + // std::cout << botC12E2Index.size() << std::endl; + } - std::cout << topC12E2Index.size() << std::endl; - std::cout << botC12E2Index.size() << std::endl; - - std::cout << topC12E2MIndex.size() << std::endl; - std::cout << botC12E2MIndex.size() << std::endl; - phiStruct phitemplate; - std::cout << "check" << std::endl; - + std::cout << topC12E2Index.size() << std::endl; + std::cout << botC12E2Index.size() << std::endl; + + std::cout << topC12E2MIndex.size() << std::endl; + std::cout << botC12E2MIndex.size() << std::endl; + phiStruct phitemplate; + std::cout << "check" << std::endl; + + for (unsigned int i = 0; i <= inputTotal.size() - 1; ++i) { + + for (unsigned phiIndex = 0; phiIndex <= 100; ++phiIndex) { + phitemplate.dist = phiIndex; + phitemplate.topPhiC12E2Count = 0; + phitemplate.botPhiC12E2Count = 0; + phitemplate.topPhiC12E2MCount = 0; + phitemplate.botPhiC12E2MCount = 0; + phiCount.push_back(phitemplate); + } - for (unsigned int i = 0; i <= inputTotal.size()-1; ++i) { - - for (unsigned phiIndex = 0; phiIndex <= 100; ++phiIndex) { - phitemplate.dist = phiIndex; - phitemplate.topPhiC12E2Count = 0; - phitemplate.botPhiC12E2Count = 0; - phitemplate.topPhiC12E2MCount = 0; - phitemplate.botPhiC12E2MCount = 0; - phiCount.push_back(phitemplate); - } - - NPX = inputTotal[i][71312].x; // x coordinate of the NP - NPY = inputTotal[i][71312].y; // y coordinate of the NP - NPZ = inputTotal[i][71312].z; // z coordinate of the NP - - int topC12E2Counter, botC12E2Counter, topC12E2MCounter, botC12E2MCounter; - - for (unsigned int dist = 0; dist <= 100; ++dist) { - - for (unsigned int index = 0; index < topC12E2Index.size(); index++) { - DistVec1 = trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][topC12E2Index[index]].x, &inputTotal[i][topC12E2Index[index]].y, &inputTotal[i][topC12E2Index[index]].z); - if (DistVec1 <= dist + 0.5 && DistVec1 >= dist - 0.5) { - phiCount[dist].topPhiC12E2Count += 1; - } - } - - for (unsigned int index = 0; index < botC12E2Index.size(); index++) { - DistVec2 = trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][botC12E2Index[index]].x, &inputTotal[i][botC12E2Index[index]].y, &inputTotal[i][botC12E2Index[index]].z); - if (DistVec2 <= dist + 0.5 && DistVec2 >= dist - 0.5) { - phiCount[dist].botPhiC12E2Count += 1; - } - } - - for (unsigned int index = 0; index < topC12E2MIndex.size(); index++) { - DistVec3 = trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][topC12E2MIndex[index]].x, &inputTotal[i][topC12E2MIndex[index]].y, &inputTotal[i][topC12E2MIndex[index]].z); - if (DistVec3 <= dist + 0.5 && DistVec3 >= dist - 0.5) { - phiCount[dist].topPhiC12E2MCount += 1; - } - } - - for (unsigned int index = 0; index < botC12E2MIndex.size(); index++) { - DistVec4 = trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][botC12E2MIndex[index]].x, &inputTotal[i][botC12E2MIndex[index]].y, &inputTotal[i][botC12E2MIndex[index]].z); - if (DistVec4 <= dist + 0.5 && DistVec4 >= dist - 0.5) { - phiCount[dist].botPhiC12E2MCount += 1; - } - } - } - phiTotal.push_back(phiCount); - phiCount.clear(); + NPX = inputTotal[i][71312].x; // x coordinate of the NP + NPY = inputTotal[i][71312].y; // y coordinate of the NP + NPZ = inputTotal[i][71312].z; // z coordinate of the NP + + int topC12E2Counter, botC12E2Counter, topC12E2MCounter, botC12E2MCounter; + + for (unsigned int dist = 0; dist <= 100; ++dist) { + + for (unsigned int index = 0; index < topC12E2Index.size(); index++) { + DistVec1 = + trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][topC12E2Index[index]].x, + &inputTotal[i][topC12E2Index[index]].y, + &inputTotal[i][topC12E2Index[index]].z); + if (DistVec1 <= dist + 0.5 && DistVec1 >= dist - 0.5) { + phiCount[dist].topPhiC12E2Count += 1; + } + } + + for (unsigned int index = 0; index < botC12E2Index.size(); index++) { + DistVec2 = + trueDist(&NPX, &NPY, &NPZ, &inputTotal[i][botC12E2Index[index]].x, + &inputTotal[i][botC12E2Index[index]].y, + &inputTotal[i][botC12E2Index[index]].z); + if (DistVec2 <= dist + 0.5 && DistVec2 >= dist - 0.5) { + phiCount[dist].botPhiC12E2Count += 1; + } + } + + for (unsigned int index = 0; index < topC12E2MIndex.size(); index++) { + DistVec3 = trueDist(&NPX, &NPY, &NPZ, + &inputTotal[i][topC12E2MIndex[index]].x, + &inputTotal[i][topC12E2MIndex[index]].y, + &inputTotal[i][topC12E2MIndex[index]].z); + if (DistVec3 <= dist + 0.5 && DistVec3 >= dist - 0.5) { + phiCount[dist].topPhiC12E2MCount += 1; + } + } + + for (unsigned int index = 0; index < botC12E2MIndex.size(); index++) { + DistVec4 = trueDist(&NPX, &NPY, &NPZ, + &inputTotal[i][botC12E2MIndex[index]].x, + &inputTotal[i][botC12E2MIndex[index]].y, + &inputTotal[i][botC12E2MIndex[index]].z); + if (DistVec4 <= dist + 0.5 && DistVec4 >= dist - 0.5) { + phiCount[dist].botPhiC12E2MCount += 1; + } + } } + phiTotal.push_back(phiCount); + phiCount.clear(); + } } - void PhiPrint () { - phipm A; + void PhiPrint() { + phipm A; for (unsigned phiIndex = 0; phiIndex <= 100; ++phiIndex) { A.phim = 0.0; A.phip = 0.0; @@ -500,197 +602,271 @@ class compute { double phi1, phi2; double phip, phim; - for (unsigned int index = 0; index <= phiTotal.size()-1; index++) { - + for (unsigned int index = 0; index <= phiTotal.size() - 1; index++) { for (unsigned int index2 = 0; index2 <= phiTotal[0].size(); index2++) { - phi1 = 0.0; - phi2 = 0.0; - phip = 0.0; - phim = 0.0; - - //std::cout << index << " " << index2 << " " << phiTotal[index][index2].topPhiC12E2Count << " " << phiTotal[index][index2].botPhiC12E2Count << " " << phiTotal[index][index2].topPhiC12E2MCount << " " << phiTotal[index][index2].botPhiC12E2MCount << std::endl; - if ((double)(phiTotal[index][index2].topPhiC12E2Count + phiTotal[index][index2].topPhiC12E2MCount) == 0.0 || (double)(phiTotal[index][index2].botPhiC12E2Count + phiTotal[index][index2].botPhiC12E2MCount) == 0) { - // Do nothing - - } else { - phi1 = (double)(phiTotal[index][index2].topPhiC12E2Count - phiTotal[index][index2].topPhiC12E2MCount)/(double)(phiTotal[index][index2].topPhiC12E2Count + phiTotal[index][index2].topPhiC12E2MCount); - phi2 = (double)(phiTotal[index][index2].botPhiC12E2Count - phiTotal[index][index2].botPhiC12E2MCount)/(double)(phiTotal[index][index2].botPhiC12E2Count + phiTotal[index][index2].botPhiC12E2MCount); - } - - phip = phi2 + phi1/2.0; - phim = phi2 - phi1/2.0; - - NewNew[index2].phim += phim; - NewNew[index2].phip += phip; - NewNew[index2].phimVec.push_back(phim); - NewNew[index2].phipVec.push_back(phip); - - std::cout << index << " " << index2 << " " << phip << " " << phim << std::endl; + phi1 = 0.0; + phi2 = 0.0; + phip = 0.0; + phim = 0.0; + + // std::cout << index << " " << index2 << " " << + // phiTotal[index][index2].topPhiC12E2Count << " " << + // phiTotal[index][index2].botPhiC12E2Count << " " << + // phiTotal[index][index2].topPhiC12E2MCount << " " << + // phiTotal[index][index2].botPhiC12E2MCount << std::endl; + if ((double)(phiTotal[index][index2].topPhiC12E2Count + + phiTotal[index][index2].topPhiC12E2MCount) == 0.0 || + (double)(phiTotal[index][index2].botPhiC12E2Count + + phiTotal[index][index2].botPhiC12E2MCount) == 0) { + // Do nothing + + } else { + phi1 = (double)(phiTotal[index][index2].topPhiC12E2Count - + phiTotal[index][index2].topPhiC12E2MCount) / + (double)(phiTotal[index][index2].topPhiC12E2Count + + phiTotal[index][index2].topPhiC12E2MCount); + phi2 = (double)(phiTotal[index][index2].botPhiC12E2Count - + phiTotal[index][index2].botPhiC12E2MCount) / + (double)(phiTotal[index][index2].botPhiC12E2Count + + phiTotal[index][index2].botPhiC12E2MCount); + } + + phip = phi2 + phi1 / 2.0; + phim = phi2 - phi1 / 2.0; + + NewNew[index2].phim += phim; + NewNew[index2].phip += phip; + NewNew[index2].phimVec.push_back(phim); + NewNew[index2].phipVec.push_back(phip); + + std::cout << index << " " << index2 << " " << phip << " " << phim + << std::endl; } } } void ComputePhiStandardDev() { double sum, mean, sq_sum, stdev; - for (std::vector::iterator it = NewNew.begin(); it != NewNew.end(); it++) { - //std::cout << " " << it - NewNew.begin() << " " << (it->phim)/inputTotal.size() << " " << (it->phip)/inputTotal.size() << " " << std::endl; - //std::cout << (it->phimVec).size() << " " << (it->phipVec).size() << std::endl; - for (unsigned int i = 0; i != NewNew[it - NewNew.begin()].phimVec.size(); ++i) { - - sum = std::accumulate(NewNew[it - NewNew.begin()].phimVec.begin(), NewNew[it - NewNew.begin()].phimVec.end(),0.0); // Compute sum - - mean = sum / NewNew[it - NewNew.begin()].phimVec.size(); // Compute Mean - - sq_sum = std::inner_product(NewNew[it - NewNew.begin()].phimVec.begin(), NewNew[it - NewNew.begin()].phimVec.end(), NewNew[it - NewNew.begin()].phimVec.begin(), 0.0); // Compute square sum - - stdev = std::sqrt(sq_sum / NewNew[it - NewNew.begin()].phimVec.size() - mean * mean)/ (pow(numberofSS,0.5)); - + for (std::vector::iterator it = NewNew.begin(); it != NewNew.end(); + it++) { + // std::cout << " " << it - NewNew.begin() << " " << + // (it->phim)/inputTotal.size() << " " << (it->phip)/inputTotal.size() << + // " " << std::endl; std::cout << (it->phimVec).size() << " " << + // (it->phipVec).size() << std::endl; + for (unsigned int i = 0; i != NewNew[it - NewNew.begin()].phimVec.size(); + ++i) { + + sum = std::accumulate(NewNew[it - NewNew.begin()].phimVec.begin(), + NewNew[it - NewNew.begin()].phimVec.end(), + 0.0); // Compute sum + + mean = sum / NewNew[it - NewNew.begin()].phimVec.size(); // Compute Mean + + sq_sum = std::inner_product(NewNew[it - NewNew.begin()].phimVec.begin(), + NewNew[it - NewNew.begin()].phimVec.end(), + NewNew[it - NewNew.begin()].phimVec.begin(), + 0.0); // Compute square sum + + stdev = std::sqrt(sq_sum / NewNew[it - NewNew.begin()].phimVec.size() - + mean * mean) / + (pow(numberofSS, 0.5)); } - - std::cout << " " << it - NewNew.begin() << " " << (it->phip)/inputTotal.size() << " " << stdev << std::endl; + + std::cout << " " << it - NewNew.begin() << " " + << (it->phip) / inputTotal.size() << " " << stdev << std::endl; } } - - void ComputeOrderphobic() { // Computes the phi, or the mismatch between the bilayer leaflets around the NP - - double DIST, DIST2; // TODO - double DISTM, DISTM2; // TODO + void ComputeOrderphobic() { // Computes the phi, or the mismatch between the + // bilayer leaflets around the NP - for (unsigned int i = 0; i < inputTotal.size(); ++i) { + double DIST, DIST2; // TODO + double DISTM, DISTM2; // TODO - for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { + for (unsigned int i = 0; i < inputTotal.size(); ++i) { - for (unsigned int newindex = 0; newindex < C12E2IndexVector.size(); ++newindex) { - - DIST = trueDist(&inputTotal[i][C12E2IndexVector[index]+4].x, &inputTotal[i][C12E2IndexVector[index]+4].y, &inputTotal[i][C12E2IndexVector[index]+4].z, &inputTotal[i][C12E2IndexVector[newindex]+4].x, &inputTotal[i][C12E2IndexVector[newindex]+4].y, &inputTotal[i][C12E2IndexVector[newindex]+4].z); - C12E2sample.index = C12E2IndexVector[index]; - C12E2sample.dist = DIST; - C12E2sample.selfXcoord = inputTotal[i][C12E2IndexVector[index]+4].x; - C12E2sample.selfYcoord = inputTotal[i][C12E2IndexVector[index]+4].y; - C12E2sample.selfZcoord = inputTotal[i][C12E2IndexVector[index]+4].z; - C12E2sample.Xcoord = inputTotal[i][C12E2IndexVector[newindex]+4].x; - C12E2sample.Ycoord = inputTotal[i][C12E2IndexVector[newindex]+4].y; - C12E2sample.Ycoord = inputTotal[i][C12E2IndexVector[newindex]+4].z; - C12E2orderphobic.push_back(C12E2sample); - - } + for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { - for (unsigned int newindex = 0; newindex < C12E2IndexVector.size(); ++newindex) { + for (unsigned int newindex = 0; newindex < C12E2IndexVector.size(); + ++newindex) { - DISTM = trueDist(&inputTotal[i][C12E2IndexVector[index]+4].x, &inputTotal[i][C12E2IndexVector[index]+4].y, &inputTotal[i][C12E2IndexVector[index]+4].z, &inputTotal[i][C12E2MIndexVector[newindex]+4].x, &inputTotal[i][C12E2MIndexVector[newindex]+4].y, &inputTotal[i][C12E2MIndexVector[newindex]+4].z); + DIST = trueDist(&inputTotal[i][C12E2IndexVector[index] + 4].x, + &inputTotal[i][C12E2IndexVector[index] + 4].y, + &inputTotal[i][C12E2IndexVector[index] + 4].z, + &inputTotal[i][C12E2IndexVector[newindex] + 4].x, + &inputTotal[i][C12E2IndexVector[newindex] + 4].y, + &inputTotal[i][C12E2IndexVector[newindex] + 4].z); C12E2sample.index = C12E2IndexVector[index]; - C12E2sample.dist = DISTM; - C12E2sample.selfXcoord = inputTotal[i][C12E2IndexVector[index]+4].x; - C12E2sample.selfYcoord = inputTotal[i][C12E2IndexVector[index]+4].y; - C12E2sample.selfZcoord = inputTotal[i][C12E2IndexVector[index]+4].z; - C12E2sample.Xcoord = inputTotal[i][C12E2MIndexVector[newindex]+4].x; - C12E2sample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex]+4].y; - C12E2sample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex]+4].z; - C12E2orderphobic.push_back(C12E2sample); - } - - sort(C12E2orderphobic.begin(), C12E2orderphobic.end(), compareByIndexOPh); - C12E2orderphobic.erase(C12E2orderphobic.begin()); - C12E2orderphobic.erase(C12E2orderphobic.begin()+6, C12E2orderphobic.end()); - // C12E2orderphobic.resize(6); - C12E2orderphobicVec.push_back(C12E2orderphobic); - - - for (unsigned int newindex = 0; newindex < C12E2MIndexVector.size(); ++newindex) { - DISTM2 = trueDist(&inputTotal[i][C12E2MIndexVector[index]+4].x, &inputTotal[i][C12E2MIndexVector[index]+4].y, &inputTotal[i][C12E2MIndexVector[index]+4].z, &inputTotal[i][C12E2MIndexVector[newindex]+4].x, &inputTotal[i][C12E2MIndexVector[newindex]+4].y, &inputTotal[i][C12E2MIndexVector[newindex]+4].z); - - C12E2Msample.index = C12E2MIndexVector[index]; - C12E2Msample.dist = DISTM2; - C12E2Msample.selfXcoord = inputTotal[i][C12E2MIndexVector[index]+4].x; - C12E2Msample.selfYcoord = inputTotal[i][C12E2MIndexVector[index]+4].y; - C12E2Msample.selfZcoord = inputTotal[i][C12E2MIndexVector[index]+4].z; - C12E2Msample.Xcoord = inputTotal[i][C12E2MIndexVector[newindex]+4].x; - C12E2Msample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex]+4].y; - C12E2Msample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex]+4].z; - C12E2Morderphobic.push_back(C12E2Msample); - - } - - - for (unsigned int newindex = 0; newindex < C12E2MIndexVector.size(); ++newindex) { - DIST2 = trueDist(&inputTotal[i][C12E2MIndexVector[index]+4].x, &inputTotal[i][C12E2MIndexVector[index]+4].y, &inputTotal[i][C12E2MIndexVector[index]+4].z, &inputTotal[i][C12E2IndexVector[newindex]+4].x, &inputTotal[i][C12E2IndexVector[newindex]+4].y, &inputTotal[i][C12E2IndexVector[newindex]+4].z); - C12E2Msample.index = C12E2MIndexVector[index]; - C12E2Msample.dist = DIST2; - C12E2Msample.selfXcoord = inputTotal[i][C12E2MIndexVector[index]+4].x; - C12E2Msample.selfYcoord = inputTotal[i][C12E2MIndexVector[index]+4].y; - C12E2Msample.selfZcoord = inputTotal[i][C12E2MIndexVector[index]+4].z; - C12E2Msample.Xcoord = inputTotal[i][C12E2IndexVector[newindex]+4].x; - C12E2Msample.Ycoord = inputTotal[i][C12E2IndexVector[newindex]+4].y; - C12E2Msample.Ycoord = inputTotal[i][C12E2IndexVector[newindex]+4].z; - C12E2Morderphobic.push_back(C12E2Msample); - } - - - //std::sort(c12E2Morderphobic.begin(), c12E2Morderphobic.end(), compareByIndexOPh); - //std::sort(C12E2Morderphobic.begin(), C12E2Morderphobic.end(), std::greater()); - //C12E2Morderphobic.erase(C12E2Morderphobic.begin()); - //C12E2Morderphobic.resize(6); - - sort(C12E2Morderphobic.begin(), C12E2Morderphobic.end(), compareByIndexOPh); - C12E2Morderphobic.erase(C12E2Morderphobic.begin()); - C12E2Morderphobic.erase(C12E2Morderphobic.begin()+6, C12E2Morderphobic.end()); - - //C12E2Morderphobic.resize(6); - - C12E2MorderphobicVec.push_back(C12E2Morderphobic); - - C12E2orderphobic.clear(); - C12E2Morderphobic.clear(); - + C12E2sample.dist = DIST; + C12E2sample.selfXcoord = inputTotal[i][C12E2IndexVector[index] + 4].x; + C12E2sample.selfYcoord = inputTotal[i][C12E2IndexVector[index] + 4].y; + C12E2sample.selfZcoord = inputTotal[i][C12E2IndexVector[index] + 4].z; + C12E2sample.Xcoord = inputTotal[i][C12E2IndexVector[newindex] + 4].x; + C12E2sample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].y; + C12E2sample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].z; + C12E2orderphobic.push_back(C12E2sample); + } + + for (unsigned int newindex = 0; newindex < C12E2IndexVector.size(); + ++newindex) { + + DISTM = trueDist(&inputTotal[i][C12E2IndexVector[index] + 4].x, + &inputTotal[i][C12E2IndexVector[index] + 4].y, + &inputTotal[i][C12E2IndexVector[index] + 4].z, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].x, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].y, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].z); + + C12E2sample.index = C12E2IndexVector[index]; + C12E2sample.dist = DISTM; + C12E2sample.selfXcoord = inputTotal[i][C12E2IndexVector[index] + 4].x; + C12E2sample.selfYcoord = inputTotal[i][C12E2IndexVector[index] + 4].y; + C12E2sample.selfZcoord = inputTotal[i][C12E2IndexVector[index] + 4].z; + C12E2sample.Xcoord = inputTotal[i][C12E2MIndexVector[newindex] + 4].x; + C12E2sample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex] + 4].y; + C12E2sample.Ycoord = inputTotal[i][C12E2MIndexVector[newindex] + 4].z; + C12E2orderphobic.push_back(C12E2sample); + } + + sort(C12E2orderphobic.begin(), C12E2orderphobic.end(), + compareByIndexOPh); + C12E2orderphobic.erase(C12E2orderphobic.begin()); + C12E2orderphobic.erase(C12E2orderphobic.begin() + 6, + C12E2orderphobic.end()); + // C12E2orderphobic.resize(6); + C12E2orderphobicVec.push_back(C12E2orderphobic); + + for (unsigned int newindex = 0; newindex < C12E2MIndexVector.size(); + ++newindex) { + DISTM2 = trueDist(&inputTotal[i][C12E2MIndexVector[index] + 4].x, + &inputTotal[i][C12E2MIndexVector[index] + 4].y, + &inputTotal[i][C12E2MIndexVector[index] + 4].z, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].x, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].y, + &inputTotal[i][C12E2MIndexVector[newindex] + 4].z); + + C12E2Msample.index = C12E2MIndexVector[index]; + C12E2Msample.dist = DISTM2; + + C12E2Msample.selfXcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].x; + C12E2Msample.selfYcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].y; + C12E2Msample.selfZcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].z; + + C12E2Msample.Xcoord = + inputTotal[i][C12E2MIndexVector[newindex] + 4].x; + C12E2Msample.Ycoord = + inputTotal[i][C12E2MIndexVector[newindex] + 4].y; + C12E2Msample.Ycoord = + inputTotal[i][C12E2MIndexVector[newindex] + 4].z; + + C12E2Morderphobic.push_back(C12E2Msample); + } + + for (unsigned int newindex = 0; newindex < C12E2MIndexVector.size(); + ++newindex) { + DIST2 = trueDist(&inputTotal[i][C12E2MIndexVector[index] + 4].x, + &inputTotal[i][C12E2MIndexVector[index] + 4].y, + &inputTotal[i][C12E2MIndexVector[index] + 4].z, + &inputTotal[i][C12E2IndexVector[newindex] + 4].x, + &inputTotal[i][C12E2IndexVector[newindex] + 4].y, + &inputTotal[i][C12E2IndexVector[newindex] + 4].z); + + C12E2Msample.index = C12E2MIndexVector[index]; + + C12E2Msample.dist = DIST2; + + C12E2Msample.selfXcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].x; + C12E2Msample.selfYcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].y; + C12E2Msample.selfZcoord = + inputTotal[i][C12E2MIndexVector[index] + 4].z; + + C12E2Msample.Xcoord = inputTotal[i][C12E2IndexVector[newindex] + 4].x; + C12E2Msample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].y; + C12E2Msample.Ycoord = inputTotal[i][C12E2IndexVector[newindex] + 4].z; + + C12E2Morderphobic.push_back(C12E2Msample); + } + + // std::sort(c12E2Morderphobic.begin(), c12E2Morderphobic.end(), + // compareByIndexOPh); std::sort(C12E2Morderphobic.begin(), + // C12E2Morderphobic.end(), std::greater()); + // C12E2Morderphobic.erase(C12E2Morderphobic.begin()); + // C12E2Morderphobic.resize(6); + + sort(C12E2Morderphobic.begin(), C12E2Morderphobic.end(), + compareByIndexOPh); + C12E2Morderphobic.erase(C12E2Morderphobic.begin()); + C12E2Morderphobic.erase(C12E2Morderphobic.begin() + 6, + C12E2Morderphobic.end()); + + // C12E2Morderphobic.resize(6); + + C12E2MorderphobicVec.push_back(C12E2Morderphobic); + + C12E2orderphobic.clear(); + C12E2Morderphobic.clear(); } - - // Final push_back + + // Final push_back orderphobicC12E2.push_back(C12E2orderphobicVec); orderphobicC12E2M.push_back(C12E2MorderphobicVec); C12E2MorderphobicVec.clear(); C12E2orderphobicVec.clear(); } - } - void OrderphobicSort() { // Computes the phi, or the mismatch between the bilayer leaflets around the NP + void OrderphobicSort() { + // Computes the phi, or the mismatch between the + // bilayer leaflets around the NP + OPVal Val; std::vector tempVal; double output; - for (unsigned int i = 0; i < orderphobicC12E2.size()-1; ++i) { - for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { - output = calcAngle(&orderphobicC12E2[i][index]); - std::cout << output << " " << i << " "<< index << " " << orderphobicC12E2[i][index][0].selfXcoord << " " << orderphobicC12E2[i][index][0].selfYcoord << " " << orderphobicC12E2[i][index][0].selfZcoord << std::endl; + for (unsigned int i = 0; i < orderphobicC12E2.size() - 1; ++i) { + for (unsigned int index = 0; index < C12E2IndexVector.size(); ++index) { + output = calcAngle(&orderphobicC12E2[i][index]); + + std::cout << output << " " << i << " " << index << " " + << orderphobicC12E2[i][index][0].selfXcoord << " " + << orderphobicC12E2[i][index][0].selfYcoord << " " + << orderphobicC12E2[i][index][0].selfZcoord << std::endl; - Val.X = orderphobicC12E2[i][index][0].selfXcoord; - Val.Y = orderphobicC12E2[i][index][0].selfYcoord; - Val.Z = orderphobicC12E2[i][index][0].selfZcoord; + Val.X = orderphobicC12E2[i][index][0].selfXcoord; + Val.Y = orderphobicC12E2[i][index][0].selfYcoord; + Val.Z = orderphobicC12E2[i][index][0].selfZcoord; Val.orderphobicVal = output; - tempVal.push_back(Val); + tempVal.push_back(Val); } - - - for (unsigned int index = 0; index < C12E2MIndexVector.size(); ++index) { - output = calcAngle(&orderphobicC12E2M[i][index]); - std::cout << output << " " << i << " "<< index << " " << orderphobicC12E2M[i][index][0].selfXcoord << " " << orderphobicC12E2M[i][index][0].selfYcoord << " " << orderphobicC12E2M[i][index][0].selfZcoord << std::endl; - Val.X = orderphobicC12E2M[i][index][0].selfXcoord; - Val.Y = orderphobicC12E2M[i][index][0].selfYcoord; - Val.Z = orderphobicC12E2M[i][index][0].selfZcoord; - Val.orderphobicVal = output; - tempVal.push_back(Val); - } - - orderphobicVectorFinal.push_back(tempVal); - tempVal.clear(); + + for (unsigned int index = 0; index < C12E2MIndexVector.size(); ++index) { + output = calcAngle(&orderphobicC12E2M[i][index]); + std::cout << output << " " << i << " " << index << " " + << orderphobicC12E2M[i][index][0].selfXcoord << " " + << orderphobicC12E2M[i][index][0].selfYcoord << " " + << orderphobicC12E2M[i][index][0].selfZcoord << std::endl; + Val.X = orderphobicC12E2M[i][index][0].selfXcoord; + Val.Y = orderphobicC12E2M[i][index][0].selfYcoord; + Val.Z = orderphobicC12E2M[i][index][0].selfZcoord; + Val.orderphobicVal = output; + tempVal.push_back(Val); + } + + orderphobicVectorFinal.push_back(tempVal); + tempVal.clear(); } } - + void printop() { + OPVal PlaceHolder; double distP; for (unsigned int dist = 0; dist <= 100; ++dist) { @@ -699,58 +875,70 @@ class compute { PlaceHolder.Z = 0.0; ABC.push_back(PlaceHolder); } - + for (unsigned int i = 0; i != orderphobicVectorFinal.size(); i++) { for (unsigned int j = 0; j != orderphobicVectorFinal[i].size(); j++) { - NPX = inputTotal[i][71312].x; // x coordinate of the NP - NPY = inputTotal[i][71312].y; // y coordinate of the NP - NPZ = inputTotal[i][71312].z; // z coordinate of the NP - - //std::cout << i << " " << j << " " << orderphobicVectorFinal[i][j].X << " " << orderphobicVectorFinal[i][j].Y << " " << orderphobicVectorFinal[i][j].Z << " " << orderphobicVectorFinal[i][j].orderphobicVal << " " << NPX << " " << NPY << " " << NPZ << std::endl; - - distP = trueDist(&NPX, &NPY, &NPZ, &orderphobicVectorFinal[i][j].X, &orderphobicVectorFinal[i][j].Y, &orderphobicVectorFinal[i][j].Z); - - for (unsigned int dist = 0; dist <= 100; ++dist) { - if (distP >= dist - 0.5 && distP <= dist + 0.5) { - ABC[dist].VV.push_back(orderphobicVectorFinal[i][j].orderphobicVal); - } - } + NPX = inputTotal[i][71312].x; // x coordinate of the NP + NPY = inputTotal[i][71312].y; // y coordinate of the NP + NPZ = inputTotal[i][71312].z; // z coordinate of the NP + + // std::cout << i << " " << j << " " << orderphobicVectorFinal[i][j].X + // << " " << orderphobicVectorFinal[i][j].Y << " " << + // orderphobicVectorFinal[i][j].Z << " " << + // orderphobicVectorFinal[i][j].orderphobicVal << " " << NPX << " " << + // NPY << " " << NPZ << std::endl; + + distP = trueDist(&NPX, &NPY, &NPZ, &orderphobicVectorFinal[i][j].X, + &orderphobicVectorFinal[i][j].Y, + &orderphobicVectorFinal[i][j].Z); + + for (unsigned int dist = 0; dist <= 100; ++dist) { + if (distP >= dist - 0.5 && distP <= dist + 0.5) { + ABC[dist].VV.push_back(orderphobicVectorFinal[i][j].orderphobicVal); + } + } } } } void LargePrint() { double sum, mean; - + for (unsigned int i = 0; i != ABC.size(); i++) { - sum = std::accumulate(ABC[i].VV.begin(), ABC[i].VV.end(),0.0); // Compute sum + sum = std::accumulate(ABC[i].VV.begin(), ABC[i].VV.end(), + 0.0); // Compute sum mean = sum / ABC[i].VV.size(); // Compute Mean - std::cout << i << " " << mean << " " << std::endl; - + std::cout << i << " " << mean << " " << std::endl; } } - + void allocateCOM() { - CenterOfMass(&C12E2IndexVector, &C12E2MIndexVector, &inputTotal, &C12E2COM, &C12E2MCOM, &C12E2TotalCOMArray, &C12E2MTotalCOMArray); + CenterOfMass(&C12E2IndexVector, &C12E2MIndexVector, &inputTotal, &C12E2COM, + &C12E2MCOM, &C12E2TotalCOMArray, &C12E2MTotalCOMArray); } -private: +private: std::vector inputVectorStruct; // push back all structs - std::vector< std::vector > > orderphobicC12E2; // This needs to be a simplified - std::vector< std::vector > > orderphobicC12E2M; // This needs to be simplfied also.. - // Vectors to store trajectory values + std::vector>> + orderphobicC12E2; // This needs to be a simplified + std::vector>> + orderphobicC12E2M; // This needs to be simplfied also.. + // Vectors to store trajectory values std::vector inputVector; // push back all structs - std::vector > inputTotal; // push back vector of structs - std::vector C12E2IndexVector; // push back C12E2 bead 7 indices (headgroups) - std::vector C12E2MIndexVector; // push back all C12E2M bead 13 indices (headgroups) - // Vectors to store the Centers of Mass + std::vector> + inputTotal; // push back vector of structs + std::vector + C12E2IndexVector; // push back C12E2 bead 7 indices (headgroups) + std::vector + C12E2MIndexVector; // push back all C12E2M bead 13 indices (headgroups) + // Vectors to store the Centers of Mass std::vector C12E2COM; std::vector C12E2MCOM; std::vector phiCount; - std::vector > phiTotal; - std::vector > C12E2TotalCOMArray; - std::vector > C12E2MTotalCOMArray; + std::vector> phiTotal; + std::vector> C12E2TotalCOMArray; + std::vector> C12E2MTotalCOMArray; // top head groups C12E2 std::vector topC12E2Index; // top head groups C12E2M @@ -759,41 +947,40 @@ class compute { std::vector topC12E2MIndex; // bot head groups C12E2M std::vector botC12E2MIndex; - // TODO - std::vector NewNew; // TODO - need to rename - - FILE *ipf; /* input file */ - int atomtype; + // TODO + std::vector NewNew; // TODO - need to rename + + FILE *ipf; /* input file */ + int atomtype; int index, l, n; - int nlines = numberofatoms + 9; - double x,y,z; /*coordinates for the atoms in the box*/ + int nlines = numberofatoms + 9; + double x, y, z; /*coordinates for the atoms in the box*/ double box1; - double box2; - double NPX, NPY, NPZ; + double box2; + double NPX, NPY, NPZ; double boxlength[boxdim]; - char line[100]; + char line[100]; inputCoord inputline; double DistVec1, DistVec2, DistVec3, DistVec4; OPHstruct C12E2sample; OPHstruct C12E2Msample; - - std::vector > C12E2orderphobicVec; // TODO - std::vector > C12E2MorderphobicVec; // TODO + + std::vector> C12E2orderphobicVec; // TODO + std::vector> C12E2MorderphobicVec; // TODO std::vector C12E2orderphobic; std::vector C12E2Morderphobic; - std::vector > orderphobicVectorFinal; // TODO + std::vector> orderphobicVectorFinal; // TODO std::vector ABC; }; -class testClass { -}; +class testClass {}; // Intitiate test class compute C12E2PhiOrderphobic; - -int main (int argc, char *argv[]) { + +int main(int argc, char *argv[]) { C12E2PhiOrderphobic.storeFile(); C12E2PhiOrderphobic.sortVectors(); @@ -807,10 +994,6 @@ int main (int argc, char *argv[]) { C12E2PhiOrderphobic.OrderphobicSort(); C12E2PhiOrderphobic.printop(); C12E2PhiOrderphobic.LargePrint(); - - return 0; -} - - - + return 0; +} diff --git a/C12E2_work/README b/C12E2_work/README index 111fa3c0..a61d6f98 100644 --- a/C12E2_work/README +++ b/C12E2_work/README @@ -1,49 +1,19 @@ -/* -------------------------------------------------------------------------------- - --------------------------------------------------------------------------------- - |Calculate the Clustering of the Mixed Lipid Bilayer through a DBSCAN Algorithm | - --------------------------------------------------------------------------------- - - How does the DBSCAN algorithm work? - - -> Point P in a cluster is a 'core' point if there are a critical number of the same type of points within a distance E. - -> Point Q is a point in the cluster if there is a clear vector towards that from any of the elements of the P vector - -> All points not reachable from any other points are outliers - - Two main algorithms are required, where there are - - --------------------------------------------------------------------------------*/ - -/* - ---------------------------------- - | DBSCAN algorithm - Pseudocode:| - ---------------------------------- -*/ - -/* - DBSCAN(D, eps, MinPts) { - C = 0 - for each point P in dataset D { - if P is visited - continue next point - mark P as visited - NeighbourPts = regionQuery(P, eps) - if Sizeof(NeighbourPts) < MinPts - mark P as NOISE - else { - C = next cluster - expandCluster(P, NeighbourPts, C, eps, MinPits) - } - } - } - */ - -/* -void DBSCAN(double *, eps, Minpts); -void regionQuery(double Distcriteria, double (*COMref) (double)) { - COMref = CenterOfMass; -} -*/ - -/* Function to find the center of mass of the C12E2 */ +Some further analysis that need to be implemented: + +- Does budding only occur in mixed surfactant systems? i.e. test with a pure bilayer + +- What is the budding ‘rate’? Plot z coordinate of the NP over time. + +- Does budding correlate with de-mixing? How does the distribution of surfactants change as budding occurs? Does a certain amount of phase-separation trigger budding? What is the distribution of surfactant A and B around the NP and does this change as budding occurs? + +- If we start the membrane in a phase-separated state, does budding still occur? (Or is it the dynamics process and starting out of equilibrium that provides the driving force?) + +- We should do repeat simulations with different random starting distributions of the different surfactant types and see if we can extract the common features of budding. + +- Is budding dependent on NP size? + +- Is budding dependent on the strength of repulsion between A and B? + + \ No newline at end of file diff --git a/C12E2_work/analysis.cxx b/C12E2_work/analysis.cxx new file mode 100644 index 00000000..178d9262 --- /dev/null +++ b/C12E2_work/analysis.cxx @@ -0,0 +1,89 @@ +#include "analysis.h" +#include +#include +#include +#include +#include + +using boost::is_any_of; + +std::vector > recorded_file; + +analysis::analysis(std::string input_file) { + /* + constructor + */ + filename = input_file; +} + +analysis::~analysis() { + std::cout << "class instance destroyed" << std::endl; +} + +void analysis::read_file() { + /* + read the file of input_file + */ + std::ifstream infile(filename); + for (std::string line; getline(infile, line);) { + // ignore all lines that start with hash + + if (line[0] == '#') { + continue; + } + + else { + boost::algorithm::split(placeholder, line, is_any_of("\t "), boost::token_compress_on); + + // store mass values + + if (placeholder[0] == "mass") { + mass_store.push_back(placeholder); + } + + // store pair coefficient values + + else if (placeholder[0] == "pair_coeff") { + pair_coeff_store.push_back(placeholder); + } + + // store bond coefficient values + + else if (placeholder[0] == "bond_coeff") { + bond_coeff_store.push_back(placeholder); + } + + // store angle coefficient values + + else if (placeholder[0] == "angle_coeff") { + angle_coeff_store.push_back(placeholder); + } + placeholder.clear(); + } + } +} + +void analysis::generate_parameter() { + /* + produce a working parameter file that + has been modified with the correct parameters + + the ones we have the keep concious about at the moment are: + + Mimic - Real + ------------ + 13-7 + 12-6 + 9-3 + 10-4 + TODO - need to change the conditional input so that the parameters can easily be modified + */ + + for (row = pair_coeff_store.begin(); row != pair_coeff_store.end(); row++) { + if (row->at(1) == "7" && row->at(2) == "7") { + std::cout << " " << row->at(0) << " " << row->at(1) << " " << row->at(2) + << " " << row->at(3) << " " << row->at(4) << " " << row->at(5) << " " << row->at(6) << std::endl; + } + } +} + diff --git a/C12E2_work/analysis.cxx~ b/C12E2_work/analysis.cxx~ new file mode 100644 index 00000000..178d9262 --- /dev/null +++ b/C12E2_work/analysis.cxx~ @@ -0,0 +1,89 @@ +#include "analysis.h" +#include +#include +#include +#include +#include + +using boost::is_any_of; + +std::vector > recorded_file; + +analysis::analysis(std::string input_file) { + /* + constructor + */ + filename = input_file; +} + +analysis::~analysis() { + std::cout << "class instance destroyed" << std::endl; +} + +void analysis::read_file() { + /* + read the file of input_file + */ + std::ifstream infile(filename); + for (std::string line; getline(infile, line);) { + // ignore all lines that start with hash + + if (line[0] == '#') { + continue; + } + + else { + boost::algorithm::split(placeholder, line, is_any_of("\t "), boost::token_compress_on); + + // store mass values + + if (placeholder[0] == "mass") { + mass_store.push_back(placeholder); + } + + // store pair coefficient values + + else if (placeholder[0] == "pair_coeff") { + pair_coeff_store.push_back(placeholder); + } + + // store bond coefficient values + + else if (placeholder[0] == "bond_coeff") { + bond_coeff_store.push_back(placeholder); + } + + // store angle coefficient values + + else if (placeholder[0] == "angle_coeff") { + angle_coeff_store.push_back(placeholder); + } + placeholder.clear(); + } + } +} + +void analysis::generate_parameter() { + /* + produce a working parameter file that + has been modified with the correct parameters + + the ones we have the keep concious about at the moment are: + + Mimic - Real + ------------ + 13-7 + 12-6 + 9-3 + 10-4 + TODO - need to change the conditional input so that the parameters can easily be modified + */ + + for (row = pair_coeff_store.begin(); row != pair_coeff_store.end(); row++) { + if (row->at(1) == "7" && row->at(2) == "7") { + std::cout << " " << row->at(0) << " " << row->at(1) << " " << row->at(2) + << " " << row->at(3) << " " << row->at(4) << " " << row->at(5) << " " << row->at(6) << std::endl; + } + } +} + diff --git a/C12E2_work/analysis.h b/C12E2_work/analysis.h new file mode 100644 index 00000000..86e45b97 --- /dev/null +++ b/C12E2_work/analysis.h @@ -0,0 +1,36 @@ +#ifndef __analysis__ +#define __analysis__ + +#include +#include +#include +#include + +// boost libraries +#include +#include + +class analysis { +public: + analysis(std::string); // constructor + void read_file(); + void generate_parameter(); + void store_scaling_parameters(); + void test_parameters(); + ~analysis(); // destructor +private: + // private variables - mostly vectors + std::string filename; + std::vector > mass_store; + std::vector > pair_coeff_store; + std::vector > bond_coeff_store; + std::vector > angle_coeff_store; + std::vector placeholder; + std::vector > placeholder_II; + // iterators + std::vector::iterator i; + std::vector::iterator stringiter; + std::vector>::iterator row; +}; + +#endif diff --git a/C12E2_work/analysis_pteros.cxx b/C12E2_work/analysis_pteros.cxx new file mode 100644 index 00000000..5ed4f84c --- /dev/null +++ b/C12E2_work/analysis_pteros.cxx @@ -0,0 +1,302 @@ +#include +#include +#include +#include +#include +// dictionary equivalent +#include // library for complex numbers +#include + +// pteros library and custom class import +#include "analysis_pteros.h" +#include +#include + +// filling in functions +bool compare_by_distance(const oph_struct &a, const oph_struct &b) { + /* + */ + return a.dist < b.dist; +} + +// filling in classes + +pteros_traj_analysis::pteros_traj_analysis(std::string input_data_gro, + std::string input_data_xtc) { + /* + read the gro and xtc files + */ + + sample_system.load(input_data_gro); + sample_system.load(input_data_xtc); + std::cout << sample_system.num_frames() << std::endl; +}; + +void pteros_traj_analysis::compute_order_parameter() { + // identify the nanoparticle + std::string id = "index "; + std::string id2 = "index "; + + pteros::Selection sel0; + pteros::Selection sel1; + pteros::Selection nanoparticle = + sample_system("name 14"); // nanoparticle selection + + // loop over each residue, compute the normal and compute + // the general representative vector and compute the angle between + for (auto &it : resid_indices) { + + sel0 = sample_system(id.append(std::to_string(it.second[0]))); + sel1 = sample_system(id2.append(std::to_string(it.second[6]))); + + // for (auto it_new = begin(it.second); it_new != end (it.second); ++it_new) + // { + // std::cout << it.first << " " << *it_new << std::endl; + // } + + // for (auto& it : ) { + + std::cout << it.first << " " << it.second[1] << " " << sel0[0].name() + << " " // sel0.center << " " + << sel1[0].name() << " " + << std::endl; // sel1.center() << std::endl; + + id = "index "; + id2 = "index "; + } + + // resid_indices; +}; + +void pteros_traj_analysis::compute_budding_rate() { + /* + compute the budding rate by comparing the z coordinate absolute + distance between the nanoparticle center from the central z coordinate + of the 'lower' bilayer + */ + std::string resname = "resid "; + pteros::Selection sel0; + pteros::Selection nanoparticle = + sample_system("name 14"); // nanoparticle selection + pteros::Selection all_surfactants( + sample_system, "name 7 6 3 4 13 12 10 9"); // surfactant selection + + // as surfactant_ids is a set, it will only store the unique + // residue ids. Any duplicated id will not be stored at all and ignored + + for (pteros::Selection::iterator it = all_surfactants.begin(); + it != all_surfactants.end(); it++) { + surfactant_ids.insert(it->resid()); + } + + // explain!! + for (set_itr = surfactant_ids.begin(); set_itr != surfactant_ids.end(); + set_itr++) { + + resname.append(std::to_string(*set_itr)); + sel0 = sample_system(resname); + resid_indices[*set_itr] = std::vector(); + + for (pteros::Selection::iterator it = sel0.begin(); it != sel0.end(); + it++) { + // create vector entry for dictionary + resid_indices[*set_itr].push_back(it->index()); + } + resname = "resid "; // reset resname string + } +} + +// void pteros_traj_analysis::compute_composition_around_np() { +// // identify the residues around the nanoparticle +// +// std::string resid = "resid "; +// std::string not_resid = "not resid "; +// +// pteros::Selection sel0; +// pteros::Selection sel1; +// +// //pteros::Selection composition_np( +// // sample_system, "within 3.0 of (name 14) and name 7 6 3 4 13 12 10 +// 9"); +// // iterator-based loop: +// +// for (set_itr = resid.begin(); set_itr != resid.end(); set_itr++) { +// not_resid.append(std::to_string(*set_itr)); +// sel0 = sample_system(not_resid); +// +// //for (pteros::Selection::iterator it = sel0.begin(); +// // it !=sel0.end(); it++) { +// +// //} +// } +// +// }; + +void pteros_traj_analysis::compute_grid(int grid_unit) { + /* + create a grid of grid_unit * grid_unit * grid_unit + */ + pteros::Grid g(grid_unit, grid_unit, grid_unit); + g.populate_periodic(sample_system.select_all()); + + for (int i = 0; i < grid_unit; ++i) + for (int j = 0; j < grid_unit; ++j) + for (int k = 0; k < grid_unit; ++k) + std::cout << g.cell(i, j, k).size() << std::endl; +}; + +void pteros_traj_analysis::compute_orderphobic_parameter() { + /* + take data from orderphobic_storage and compute the + orderphobic parameter for that particular surfactant + + phi_l = |(1/6) \sum_{allpairs} exp(6 * i * theta_{ij})|^2 + + where theta is the angle between an arbitrary axis and a vector + connecting tail-end particle i to tail-end particle j + + */ + + std::complex inital_complex(0.0, 0.0); + // initalize a complex fraction 1/6 with a 0 in the complex plane + std::complex inital_sixth((double)1 / 6, 0.0); + + for (int i = 0; i <= orderphobic_storage.size(); i++) { + for (std::vector::iterator it = + orderphobic_storage.at(i).begin(); + it != orderphobic_storage.at(i).end(); ++it) { + + + } + } + +}; + +void pteros_traj_analysis::collect_orderphobic_parameter_vector() { + /* + compute the orderphobic parameter based on an arbitrary + vector pointing out and + */ + double u, v; + double dist; + // assign arbitrary vector to compare with + double reference_vector[2]; + // initialize a complex number with 0s + // in both the real and complex plane + std::complex inital_complex(0.0, 0.0); + // initalize a complex fraction 1/6 with a 0 in the complex plane + std::complex inital_sixth((double)1 / 6, 0.0); + // make sure to exclude waters + + pteros::Selection new_system; + pteros::Selection new_system_II; + pteros::Selection adhoc_system; + + reference_vector[0] = 10.0; + reference_vector[1] = 10.0; + u = pow((pow(reference_vector[0], 2) + pow(reference_vector[1], 2)), 0.5); + + // sample_system("name 14"); // nanoparticle selection + std::string adhoc_resname = "not name 1 and resid "; + oph_struct orderphobic_data; + std::vector single_storage; + + std::string resid_1 = "resid "; + std::string resid_2 = "resid "; + + // where to store orderphobic data + int count = 0; + for (auto &it : resid_indices) + + if (count == 50) { + break; + } else { + resid_1.append(std::to_string(it.first)); // get residue id + new_system = sample_system(resid_1); + + for (auto &it_2 : resid_indices) { + + if (it.first != it_2.first) { + + resid_2.append(std::to_string(it_2.first)); // get residue id + new_system_II = sample_system(resid_2); + + dist = this->cartesian_distance( + new_system_II[4].x(), new_system_II[4].y(), new_system_II[4].z(), + new_system[4].x(), new_system[4].y(), new_system[4].z()); + + // std::cout << resid_1 << " " << resid_2 << " " << + // new_system_II[4].name() << " " << new_system[4].name() << " " << + // dist << std::endl; + + // --- storing the orderphobic data --- + orderphobic_data.resid_self = it.first; + orderphobic_data.resid_other = it_2.first; + orderphobic_data.dist = dist; + // store cartesian coordinates of the self + orderphobic_data.self_x_coord = new_system_II[4].x(); + orderphobic_data.self_y_coord = new_system_II[4].y(); + orderphobic_data.self_z_coord = new_system_II[4].z(); + // store cartesian coordinates of the cooordinates away + orderphobic_data.x_coord = new_system[4].x(); + orderphobic_data.y_coord = new_system[4].y(); + orderphobic_data.z_coord = new_system[4].z(); + orderphobic_data.orderphobic_value_self = + 0.0; // placeholder for the orderphobic value + // append entry + single_storage.push_back(orderphobic_data); + // reset strings + // std::cout << newcount << std::endl; + resid_2 = "resid "; + } + } + + resid_1 = "resid "; + // sort entry and then push_back + std::sort(single_storage.begin(), single_storage.end(), + compare_by_distance); + + orderphobic_storage.push_back(single_storage); + // clear all entries for next iteration + single_storage.clear(); + count += 1; + } +}; + +void pteros_traj_analysis::print_structs(int limit = 5) { + /* + print out vector of vector entries for the orderphobic struct inputs + */ + for (int i = 0; i <= orderphobic_storage.size(); i++) { + for (std::vector::iterator it = + orderphobic_storage.at(i).begin(); + it != orderphobic_storage.at(i).end(); ++it) { + std::cout << i << " " << it->resid_self << " " << it->resid_other << " " + << it->dist << std::endl; + // break if conditions + if (it - orderphobic_storage.at(i).begin() == limit) { + break; + } + } + } +}; + +double pteros_traj_analysis::cartesian_distance(float &com_1x, float &com_1y, + float &com_1z, float &com_2x, + float &com_2y, float &com_2z) { + /* + compute absolute distance between 3d cartesian points + */ + + double cdist = pow((pow(com_1x - com_2x, 2) + pow(com_1y - com_2y, 2) + + pow(com_1z - com_2z, 2)), + 0.5); + return cdist; +}; + +pteros_traj_analysis::~pteros_traj_analysis() { + /* + class destructor + */ + std::cout << "destructor called" << std::endl; +}; diff --git a/C12E2_work/analysis_pteros.cxx~ b/C12E2_work/analysis_pteros.cxx~ new file mode 100644 index 00000000..12aae42a --- /dev/null +++ b/C12E2_work/analysis_pteros.cxx~ @@ -0,0 +1,272 @@ +#include +#include +#include +#include +#include +// dictionary equivalent +#include // library for complex numbers +#include + +// pteros library and custom class import +#include "analysis_pteros.h" +#include +#include + +// filling in functions +bool compare_by_distance(const oph_struct &a, const oph_struct &b) { + /* + */ + return a.dist < b.dist; +} + +// filling in classes + +pteros_traj_analysis::pteros_traj_analysis(std::string input_data_gro, + std::string input_data_xtc) { + /* + read the gro and xtc files + */ + + sample_system.load(input_data_gro); + sample_system.load(input_data_xtc); + std::cout << sample_system.num_frames() << std::endl; +}; + +void pteros_traj_analysis::compute_order_parameter() { + // identify the nanoparticle + std::string id = "index "; + std::string id2 = "index "; + + pteros::Selection sel0; + pteros::Selection sel1; + pteros::Selection nanoparticle = + sample_system("name 14"); // nanoparticle selection + + // loop over each residue, compute the normal and compute + // the general representative vector and compute the angle between + for (auto &it : resid_indices) { + + sel0 = sample_system(id.append(std::to_string(it.second[0]))); + sel1 = sample_system(id2.append(std::to_string(it.second[6]))); + + // for (auto it_new = begin(it.second); it_new != end (it.second); ++it_new) + // { + // std::cout << it.first << " " << *it_new << std::endl; + // } + + // for (auto& it : ) { + + std::cout << it.first << " " << it.second[1] << " " << sel0[0].name() + << " " // sel0.center << " " + << sel1[0].name() << " " + << std::endl; // sel1.center() << std::endl; + + id = "index "; + id2 = "index "; + } + + // resid_indices; +}; + +void pteros_traj_analysis::compute_budding_rate() { + /* + compute the budding rate by comparing the z coordinate absolute + distance between the nanoparticle center from the central z coordinate + of the 'lower' bilayer + */ + std::string resname = "resid "; + pteros::Selection sel0; + pteros::Selection nanoparticle = + sample_system("name 14"); // nanoparticle selection + pteros::Selection all_surfactants( + sample_system, "name 7 6 3 4 13 12 10 9"); // surfactant selection + + // as surfactant_ids is a set, it will only store the unique + // residue ids. Any duplicated id will not be stored at all and ignored + + for (pteros::Selection::iterator it = all_surfactants.begin(); + it != all_surfactants.end(); it++) { + surfactant_ids.insert(it->resid()); + } + + // explain!! + for (set_itr = surfactant_ids.begin(); set_itr != surfactant_ids.end(); + set_itr++) { + + resname.append(std::to_string(*set_itr)); + sel0 = sample_system(resname); + resid_indices[*set_itr] = std::vector(); + + for (pteros::Selection::iterator it = sel0.begin(); it != sel0.end(); + it++) { + // create vector entry for dictionary + resid_indices[*set_itr].push_back(it->index()); + } + resname = "resid "; // reset resname string + } +} + +// void pteros_traj_analysis::compute_composition_around_np() { +// // identify the residues around the nanoparticle +// +// std::string resid = "resid "; +// std::string not_resid = "not resid "; +// +// pteros::Selection sel0; +// pteros::Selection sel1; +// +// //pteros::Selection composition_np( +// // sample_system, "within 3.0 of (name 14) and name 7 6 3 4 13 12 10 +// 9"); +// // iterator-based loop: +// +// for (set_itr = resid.begin(); set_itr != resid.end(); set_itr++) { +// not_resid.append(std::to_string(*set_itr)); +// sel0 = sample_system(not_resid); +// +// //for (pteros::Selection::iterator it = sel0.begin(); +// // it !=sel0.end(); it++) { +// +// //} +// } +// +// }; + +void pteros_traj_analysis::compute_grid(int grid_unit) { + /* + create a grid of grid_unit * grid_unit * grid_unit + */ + pteros::Grid g(grid_unit, grid_unit, grid_unit); + g.populate_periodic(sample_system.select_all()); + + for (int i = 0; i < grid_unit; ++i) + for (int j = 0; j < grid_unit; ++j) + for (int k = 0; k < grid_unit; ++k) + std::cout << g.cell(i, j, k).size() << std::endl; +}; + +void pteros_traj_analysis::compute_orderphobic_parameter() { + /* + compute the orderphobic parameter based on an arbitrary + vector pointing out and + */ + double u, v; + double dist; + // assign arbitrary vector to compare with + double reference_vector[2]; + // initialize a complex number with 0s + // in both the real and complex plane + std::complex inital_complex(0.0, 0.0); + // initalize a complex fraction 1/6 with a 0 in the complex plane + std::complex inital_sixth((double)1 / 6, 0.0); + // make sure to exclude waters + + pteros::Selection new_system; + pteros::Selection new_system_II; + pteros::Selection adhoc_system; + + reference_vector[0] = 10.0; + reference_vector[1] = 10.0; + u = pow((pow(reference_vector[0], 2) + pow(reference_vector[1], 2)), 0.5); + + // sample_system("name 14"); // nanoparticle selection + std::string adhoc_resname = "not name 1 and resid "; + oph_struct orderphobic_data; + std::vector single_storage; + + std::string resid_1 = "resid "; + std::string resid_2 = "resid "; + + // where to store orderphobic data + int count = 0; + for (auto &it : resid_indices) + + if (count == 50) { + break; + } else { + resid_1.append(std::to_string(it.first)); // get residue id + new_system = sample_system(resid_1); + + for (auto &it_2 : resid_indices) { + + if (it.first != it_2.first) { + + resid_2.append(std::to_string(it_2.first)); // get residue id + new_system_II = sample_system(resid_2); + + dist = this->cartesian_distance( + new_system_II[4].x(), new_system_II[4].y(), new_system_II[4].z(), + new_system[4].x(), new_system[4].y(), + new_system[4].z()); + + //std::cout << resid_1 << " " << resid_2 << " " << new_system_II[4].name() << " " << new_system[4].name() << " " << dist << std::endl; + + // --- storing the orderphobic data --- + orderphobic_data.resid_self = it.first; + orderphobic_data.resid_other = it_2.first; + orderphobic_data.dist = dist; + // store cartesian coordinates of the self + orderphobic_data.self_x_coord = new_system_II[4].x(); + orderphobic_data.self_y_coord = new_system_II[4].y(); + orderphobic_data.self_z_coord = new_system_II[4].z(); + // store cartesian coordinates of the cooordinates away + orderphobic_data.x_coord = new_system[4].x(); + orderphobic_data.y_coord = new_system[4].y(); + orderphobic_data.z_coord = new_system[4].z(); + // append entry + single_storage.push_back(orderphobic_data); + // reset strings + //std::cout << newcount << std::endl; + resid_2 = "resid "; + } + } + + resid_1 = "resid "; + // sort entry and then push_back + std::sort(single_storage.begin(), single_storage.end(), + compare_by_distance); + + orderphobic_storage.push_back(single_storage); + // clear all entries for next iteration + single_storage.clear(); + count += 1; + } +}; + +void pteros_traj_analysis::print_structs(int limit = 5) { + /* + print out vector of vector entries for the orderphobic struct inputs + */ + for (int i = 0; i <= orderphobic_storage.size(); i++) { + for (std::vector::iterator it = orderphobic_storage.at(i).begin() ; it != orderphobic_storage.at(i).end(); ++it) { + std::cout << i << " " << it->resid_self + << " " << it->resid_other << " " + << it->dist << std::endl; + // break if conditions + if (it - orderphobic_storage.at(i).begin() == limit) { + break; + } + + } + } +}; + +double pteros_traj_analysis::cartesian_distance(float &com_1x, float &com_1y, + float &com_1z, float &com_2x, + float &com_2y, float &com_2z) { + /* + compute absolute distance between 3d cartesian points + */ + + double cdist = pow((pow(com_1x - com_2x, 2) + pow(com_1y - com_2y, 2) + + pow(com_1z - com_2z, 2)), + 0.5); + return cdist; +}; + +pteros_traj_analysis::~pteros_traj_analysis() { + /* + class destructor + */ + std::cout << "destructor called" << std::endl; +}; diff --git a/C12E2_work/analysis_pteros.h b/C12E2_work/analysis_pteros.h new file mode 100644 index 00000000..928c8fb7 --- /dev/null +++ b/C12E2_work/analysis_pteros.h @@ -0,0 +1,76 @@ +#ifndef __analysis_pteros__ +#define __analysis_pteros__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include // library for complex numbers + + +typedef struct { // Used to input the center of masses for each lipid + double dist; + int top_phi_C12E2_count; + int top_phi_C12E2M_count; + int bot_phi_C12E2_count; + int bot_phi_C12E2M_count; +} phi_struct; + +typedef struct { // Used to identify the group and distance to compute the + // orderphobic effect + float resid_self; + float resid_other; + float dist; + float self_x_coord; + float self_y_coord; + float self_z_coord; + float x_coord; + float y_coord; + float z_coord; + float orderphobic_value_self = 0.0; +} oph_struct; + + +// functions to work with +bool compare_by_distance(const oph_struct&, const oph_struct&); + +// classes to work with +class pteros_traj_analysis { +public: + pteros_traj_analysis(std::string, std::string); // constructor + ~pteros_traj_analysis(); // destructor + // void methodologies + void compute_order_parameter(); + void compute_budding_rate(); + //void compute_composition_around_np(); + void compute_grid(int); + void compute_resid(int); // new + + // orderphobic parameter collection and computation + void collect_orderphobic_parameter_vector(); + void compute_orderphobic_parameter(); + // protected methods are accesible in the class that defines them and in classes + // that inherit from that class. + double cartesian_distance(float&, float&, float&, float&, + float&, float&); + + // print functions + void print_structs(int); + + +private: // private members are only accessible within the class defining them + + std::vector > orderphobic_storage; + std::set resid; + std::set surfactant_ids; + std::set::iterator set_itr; + // unordered map - the equivalent of a dictionary + std::unordered_map > resid_indices; + pteros::System sample_system; +}; + +#endif diff --git a/C12E2_work/analysis_pteros.h~ b/C12E2_work/analysis_pteros.h~ new file mode 100644 index 00000000..5e925776 --- /dev/null +++ b/C12E2_work/analysis_pteros.h~ @@ -0,0 +1,69 @@ +#ifndef __analysis_pteros__ +#define __analysis_pteros__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include // library for complex numbers + + +typedef struct { // Used to input the center of masses for each lipid + double dist; + int top_phi_C12E2_count; + int top_phi_C12E2M_count; + int bot_phi_C12E2_count; + int bot_phi_C12E2M_count; +} phi_struct; + +typedef struct { // Used to identify the group and distance to compute the + // orderphobic effect + float resid_self; + float resid_other; + float dist; + float self_x_coord; + float self_y_coord; + float self_z_coord; + float x_coord; + float y_coord; + float z_coord; +} oph_struct; + +bool compare_by_distance(const oph_struct&, const oph_struct&); + +class pteros_traj_analysis { +public: + pteros_traj_analysis(std::string, std::string); // constructor + ~pteros_traj_analysis(); // destructor + // void methodologies + void compute_order_parameter(); + void compute_budding_rate(); + //void compute_composition_around_np(); + void compute_grid(int); + void compute_resid(int); // new + void compute_orderphobic_parameter(); + // protected methods are accesible in the class that defines them and in classes + // that inherit from that class. + double cartesian_distance(float&, float&, float&, float&, + float&, float&); + + // print functions + void print_structs(int); + + +private: // private members are only accessible within the class defining them + + std::vector > orderphobic_storage; + std::set resid; + std::set surfactant_ids; + std::set::iterator set_itr; + // unordered map - the equivalent of a dictionary + std::unordered_map > resid_indices; + pteros::System sample_system; +}; + +#endif diff --git a/C12E2_work/main.cxx b/C12E2_work/main.cxx new file mode 100644 index 00000000..7ab3ec31 --- /dev/null +++ b/C12E2_work/main.cxx @@ -0,0 +1,33 @@ +#include +#include +#include +#include + +// boost libraries +#include +#include + +// custom library +#include "analysis.h" +#include "analysis_pteros.h" + +int main(void) { + // edit forcefield files + //analysis initial("/home/sang/Desktop/git/Research/C12E2_work/parameters/parm.cg_cmm_fake"); + //initial.read_file(); + //initial.generate_parameter(); + + pteros_traj_analysis pteros_data( + "/home/sang/Desktop/git/Research/C12E2_work/trajectories/dump.gro", + "/home/sang/Desktop/git/Research/C12E2_work/trajectories/dump.xtc"); + + //pteros_data.compute_order_parameter(); + //pteros_data.compute_composition_around_np(); + //pteros_data.compute_grid(100); + //pteros_data.compute_order_parameter(); + pteros_data.compute_budding_rate(); + pteros_data.collect_orderphobic_parameter_vector(); + pteros_data.print_structs(50); + + return 0; +} diff --git a/C12E2_work/main.cxx~ b/C12E2_work/main.cxx~ new file mode 100644 index 00000000..04cd1b80 --- /dev/null +++ b/C12E2_work/main.cxx~ @@ -0,0 +1,33 @@ +#include +#include +#include +#include + +// boost libraries +#include +#include + +// custom library +#include "analysis.h" +#include "analysis_pteros.h" + +int main(void) { + // edit forcefield files + //analysis initial("/home/sang/Desktop/git/Research/C12E2_work/parameters/parm.cg_cmm_fake"); + //initial.read_file(); + //initial.generate_parameter(); + + pteros_traj_analysis pteros_data( + "/home/sang/Desktop/git/Research/C12E2_work/trajectories/dump.gro", + "/home/sang/Desktop/git/Research/C12E2_work/trajectories/dump.xtc"); + + //pteros_data.compute_order_parameter(); + //pteros_data.compute_composition_around_np(); + //pteros_data.compute_grid(100); + pteros_data.compute_budding_rate(); + //pteros_data.compute_order_parameter(); + pteros_data.compute_orderphobic_parameter(); + pteros_data.print_structs(50); + + return 0; +} diff --git a/C12E2_work/new b/C12E2_work/new new file mode 100755 index 00000000..b5df9c70 Binary files /dev/null and b/C12E2_work/new differ diff --git a/C12E2_work/new.c b/C12E2_work/new.c new file mode 100644 index 00000000..e69de29b diff --git a/C12E2_work/old.cxx b/C12E2_work/old.cxx deleted file mode 100644 index bba3a849..00000000 --- a/C12E2_work/old.cxx +++ /dev/null @@ -1,189 +0,0 @@ - -/* - - void AllocationVec() { - The way that the C12E2 and C12E2-M batches are allocated in the simulation follows - the fact that this was directly replicated in LAMMPS; hence, there isn't a single - block of list, but rather 4 batches of C12E2 and C12E2-M portions that we need to - allocate - - Each batch has 250 molecules, which is why we loop from 0 to 249 - - for (int i = 0; i <= 249; i++) { - - // C12E2 batches - - int batch1 = 0; - int batch2 = 3500; - int batch3 = 7000; - int batch4 = 10500; - - // C12E2-M batches - - int batch1M = 1750; - int batch2M = 5250; - int batch3M = 8750; - int batch4M = 12250; - - // Allocate index - - - - A7 = 7*(i); - A6_1= 7*(i) + 1; - A6_2= 7*(i) + 2; - A3_1 = 7*(i) + 3; - A3_2 = 7*(i) + 4; - A3_3 = 7*(i) + 5; - A4 = 7*(i) + 6; - - // Second batch - - A7_2 = 7*(i) + 3501; - A6_1_2 = 7*(i) + 3502; - A6_2_2 = 7*(i) + 3503; - A3_1_2 = 7*(i) + 3504; - A3_2_2 = 7*(i) + 3505; - A3_3_2 = 7*(i) + 3506; - A4_2 = 7*(i) + 3507; - - // Third batch - - A7_3 = 7*(i) + 7001; - A6_1_3 = 7*(i) + 7002; - A6_2_3 = 7*(i) + 7003; - A3_1_3 = 7*(i) + 7004; - A3_2_3 = 7*(i) + 7005; - A3_3_3 = 7*(i) + 7006; - A4_3 = 7*(i) + 7007; - - // Fourth batch - - A7_4 = 7*(i) + 10501; - A6_1_4 = 7*(i) + 10502; - A6_2_4 = 7*(i) + 10503; - A3_1_4 = 7*(i) + 10504; - A3_2_4 = 7*(i) + 10505; - A3_3_4 = 7*(i) + 10506; - A4_4 = 7*(i) + 10507; - - // --- C12E2-M --- // - - - A13 = 7*(i) + 1751; - A12_1= 7*(i)+ 1752; - A12_2= 7*(i)+ 1753; - A9_1 = 7*(i)+ 1754; - A9_2 = 7*(i)+ 1755; - A9_3 = 7*(i)+ 1756; - A10 = 7*(i)+ 1757; - - /* Mimic atoms - Second batch */ - - A13_2 = 7*(i) + 5251; - A12_1_2 = 7*(i)+ 5252; - A12_2_2 = 7*(i)+ 5253; - A9_1_2 = 7*(i)+ 5254; - A9_2_2 = 7*(i)+ 5255; - A9_3_2 = 7*(i)+ 5256; - A10_2 = 7*(i)+ 5257; - -// Mimic atoms - Third batch - - A13_3 = 7*(i) + 8751; - A12_1_3 = 7*(i)+ 8752; - A12_2_3 = 7*(i)+ 8753; - A9_1_3 = 7*(i)+ 8754; - A9_2_3 = 7*(i)+ 8755; - A9_3_3 = 7*(i)+ 8756; - A10_3 = 7*(i)+ 8757; - -// Mimic atoms - Fourth batch */ - - A13_4 = 7*(i) + 12251; - A12_1_4 = 7*(i)+ 12252; - A12_2_4 = 7*(i)+ 12253; - A9_1_4 = 7*(i)+ 12254; - A9_2_4 = 7*(i)+ 12255; - A9_3_4 = 7*(i)+ 12256; - A10_4 = 7*(i)+ 12257; - - - - // First batch - C12E2_struct[i].index[0] = A7; - C12E2_struct[i].index[1] = A6_1; - C12E2_struct[i].index[2] = A6_2; - C12E2_struct[i].index[3] = A3_1; - C12E2_struct[i].index[4] = A3_2; - C12E2_struct[i].index[5] = A3_3; - C12E2_struct[i].index[6] = A4; - - // Second batch - C12E2_struct[i+250].index[0] = A7_2; - C12E2_struct[i+250].index[1] = A6_1_2; - C12E2_struct[i+250].index[2] = A6_2_2; - C12E2_struct[i+250].index[3] = A3_1_2; - C12E2_struct[i+250].index[4] = A3_2_2; - C12E2_struct[i+250].index[5] = A3_3_2; - C12E2_struct[i+250].index[6] = A4_2; - - // Third batch - C12E2_struct[i+500].index[0] = A7_3; - C12E2_struct[i+500].index[1] = A6_1_3; - C12E2_struct[i+500].index[2] = A6_2_3; - C12E2_struct[i+500].index[3] = A3_1_3; - C12E2_struct[i+500].index[4] = A3_2_3; - C12E2_struct[i+500].index[5] = A3_3_3; - C12E2_struct[i+500].index[6] = A4_3; - - // Fourth batch - C12E2_struct[i+750].index[0] = A7_4; - C12E2_struct[i+750].index[1] = A6_1_4; - C12E2_struct[i+750].index[2] = A6_2_4; - C12E2_struct[i+750].index[3] = A3_1_4; - C12E2_struct[i+750].index[4] = A3_2_4; - C12E2_struct[i+750].index[5] = A3_3_4; - C12E2_struct[i+750].index[6] = A4_4; - - // Mimic assignment - - - // First batch - C12E2M_struct[i].index[0] = A13; - C12E2M_struct[i].index[1] = A12_1; - C12E2M_struct[i].index[2] = A12_2; - C12E2M_struct[i].index[3] = A9_1; - C12E2M_struct[i].index[4] = A9_2; - C12E2M_struct[i].index[5] = A9_3; - C12E2M_struct[i].index[6] = A10; - - // Second batch - C12E2M_struct[i+250].index[0] = A13_2; - C12E2M_struct[i+250].index[1] = A12_1_2; - C12E2M_struct[i+250].index[2] = A12_2_2; - C12E2M_struct[i+250].index[3] = A9_1_2; - C12E2M_struct[i+250].index[4] = A9_2_2; - C12E2M_struct[i+250].index[5] = A9_3_2; - C12E2M_struct[i+250].index[6] = A10_2; - - // Third batch - C12E2M_struct[i+500].index[0] = A13_3; - C12E2M_struct[i+500].index[1] = A12_1_3; - C12E2M_struct[i+500].index[2] = A12_2_3; - C12E2M_struct[i+500].index[3] = A9_1_3; - C12E2M_struct[i+500].index[4] = A9_2_3; - C12E2M_struct[i+500].index[5] = A9_3_3; - C12E2M_struct[i+500].index[6] = A10_3; - - // Fourth batch - C12E2M_struct[i+750].index[0] = A13_4; - C12E2M_struct[i+750].index[1] = A12_1_4; - C12E2M_struct[i+750].index[2] = A12_2_4; - C12E2M_struct[i+750].index[3] = A9_1_4; - C12E2M_struct[i+750].index[4] = A9_2_4; - C12E2M_struct[i+750].index[5] = A9_3_4; - C12E2M_struct[i+750].index[6] = A10_4; - } - } -*/ diff --git a/C12E2_work/parameters/parm.cg_cmm_fake b/C12E2_work/parameters/parm.cg_cmm_fake new file mode 100755 index 00000000..f5821ead --- /dev/null +++ b/C12E2_work/parameters/parm.cg_cmm_fake @@ -0,0 +1,196 @@ +## Parameter file for CG-CMM coarse grain force field +## W Shinoda et al, Mol Sim, 33, 25 (2007) +## W Shinoda et al, Soft Matter, 4, 2454 (2008) +## +## Definition of CG particle types +## 1 (W): 3H2O +## 2 (CT): CH3-CH2-CH2 +## 3 (CM): CH2-CH2-CH2 +## 4 (CT2): CH3-CH2 +## 5 (EOT): CH3-O-CH2 +## 6 (EO): CH2-O-CH2 +## 7 (OA): CH2-OH +## 8 (CFT): CF3-CF2-CF2 \ +## 9 (CF): CF2-CF2-CF2 | NB masses haven't been changed +## 10 (CFT2): CF3-CF2 / +## 11 and 12 are the equivalents to 5 and 6 +## 11 (EOTF): CF3-O-CF2 +## 12 (EOF): CF2-O-CF2 +## 13 (OAF): CF2-OH +## +## Notes: +## W-EOT non-bonded interaction parameter missing from published force field (assume same as W-EO) +## For new headgroups, 11,12 and 13, there is a decrease of 0.4 to the epsilon value compared to the original 5,6,and 7 headgroups + +#pair_style hybrid cg/cmm 15.0 table linear 2000 +#pair_style hybrid lj/sdk 15.0 table linear 2000 + +#bond_style harmonic +#angle_style sdk +#special_bonds lj/coul 0.0 0.0 1.0 + +mass 1 54.000000 +mass 2 43.089000 +mass 3 42.081000 +mass 4 29.062000 +mass 5 45.059400 +mass 6 44.051400 +mass 7 31.034400 +mass 8 43.089000 +mass 9 42.081000 +mass 10 29.062000 +mass 11 45.059400 +mass 12 44.051400 +mass 13 31.034400 +mass 14 139.9056 + +## i j pot type epsilon sigma + +pair_coeff 1 1 lj/sdk lj12_4 0.895000 4.371000 +pair_coeff 1 2 lj/sdk lj12_4 0.360000 4.478000 +pair_coeff 1 3 lj/sdk lj12_4 0.340000 4.438500 +pair_coeff 1 4 lj/sdk lj12_4 0.290000 4.296000 +pair_coeff 1 5 lj/sdk lj12_4 0.570000 4.310000 +pair_coeff 1 6 lj/sdk lj12_4 0.570000 4.310000 +pair_coeff 1 7 lj/sdk lj12_4 0.700000 3.950000 +pair_coeff 1 8 lj/sdk lj12_4 0.360000 4.478000 +pair_coeff 1 9 lj/sdk lj12_4 0.340000 4.438500 +pair_coeff 1 10 lj/sdk lj12_4 0.290000 4.296000 +pair_coeff 1 11 lj/sdk lj12_4 0.570000 4.310000 +pair_coeff 1 12 lj/sdk lj12_4 0.570000 4.310000 +pair_coeff 1 13 lj/sdk lj12_4 0.700000 3.950000 +pair_coeff 2 2 lj/sdk lj9_6 0.469000 4.585000 +pair_coeff 2 3 lj/sdk lj9_6 0.444000 4.545500 +pair_coeff 2 4 lj/sdk lj9_6 0.382500 4.403000 +pair_coeff 2 5 lj/sdk lj9_6 0.452700 4.417500 +pair_coeff 2 6 lj/sdk lj9_6 0.410000 4.340000 +pair_coeff 2 7 lj/sdk lj9_6 0.437200 4.033000 +## the epsilon values for the following three interactions should be set +## to those of the 2-2, 2-3 and 2-4 interactions, scaled by an appropriate factor +pair_coeff 2 8 lj/sdk lj9_6 0.422100 4.585000 +pair_coeff 2 9 lj/sdk lj9_6 0.399600 4.545500 +pair_coeff 2 10 lj/sdk lj9_6 0.344250 4.403000 +pair_coeff 2 11 lj/sdk lj9_6 0.452700 4.417500 +pair_coeff 2 12 lj/sdk lj9_6 0.410000 4.340000 +pair_coeff 2 13 lj/sdk lj9_6 0.437000 4.033000 +pair_coeff 3 3 lj/sdk lj9_6 0.420000 4.506000 +pair_coeff 3 4 lj/sdk lj9_6 0.362000 4.363500 +pair_coeff 3 5 lj/sdk lj9_6 0.428400 4.378000 +pair_coeff 3 6 lj/sdk lj9_6 0.377000 4.274000 +pair_coeff 3 7 lj/sdk lj9_6 0.365000 3.987000 +## the epsilon values for the following three interactions should be set +## to those of the 2-3,3-3, and 3-4 interactions, scaled by an appropriate factor +pair_coeff 3 8 lj/sdk lj9_6 0.399600 4.545500 +pair_coeff 3 9 lj/sdk lj9_6 0.390000 4.506000 ### +pair_coeff 3 10 lj/sdk lj9_6 0.332000 4.363500 ### +pair_coeff 3 11 lj/sdk lj9_6 0.428400 4.378000 +pair_coeff 3 12 lj/sdk lj9_6 0.377000 4.274000 +pair_coeff 3 13 lj/sdk lj9_6 0.365000 3.987000 +pair_coeff 4 4 lj/sdk lj9_6 0.312000 4.221000 +pair_coeff 4 5 lj/sdk lj9_6 0.369200 4.235500 +pair_coeff 4 6 lj/sdk lj9_6 0.370000 4.274000 +pair_coeff 4 7 lj/sdk lj9_6 0.380000 3.840000 +## the epsilon values for the following three interactions should be set +## to those of the 2-4,3-4, and 4-4 interactions, scaled by an appropriate factor +pair_coeff 4 8 lj/sdk lj9_6 0.344250 4.403000 +pair_coeff 4 9 lj/sdk lj9_6 0.332000 4.363500 ### +pair_coeff 4 10 lj/sdk lj9_6 0.282000 4.221000 ### +pair_coeff 4 11 lj/sdk lj9_6 0.369200 4.235500 +pair_coeff 4 12 lj/sdk lj9_6 0.370000 4.274000 +pair_coeff 4 13 lj/sdk lj9_6 0.380000 3.840000 +pair_coeff 5 5 lj/sdk lj9_6 0.437000 4.250000 +pair_coeff 5 6 lj/sdk lj9_6 0.420000 4.250000 +pair_coeff 5 7 lj/sdk lj9_6 0.443000 3.981500 +pair_coeff 5 8 lj/sdk lj9_6 0.452700 4.417500 +pair_coeff 5 9 lj/sdk lj9_6 0.428400 4.378000 +pair_coeff 5 10 lj/sdk lj9_6 0.369200 4.235500 +pair_coeff 5 11 lj/sdk lj9_6 0.397000 4.250000 +pair_coeff 5 12 lj/sdk lj9_6 0.380000 4.250000 +pair_coeff 5 13 lj/sdk lj9_6 0.403000 3.981500 +pair_coeff 6 6 lj/sdk lj9_6 0.405000 4.250000 +pair_coeff 6 7 lj/sdk lj9_6 0.440000 3.890000 +pair_coeff 6 8 lj/sdk lj9_6 0.410000 4.340000 +pair_coeff 6 9 lj/sdk lj9_6 0.377000 4.274000 +pair_coeff 6 10 lj/sdk lj9_6 0.370000 4.274000 +pair_coeff 6 11 lj/sdk lj9_6 0.380000 4.250000 +pair_coeff 6 12 lj/sdk lj9_6 0.375000 4.250000 ### +pair_coeff 6 13 lj/sdk lj9_6 0.410000 3.890000 ### +pair_coeff 7 7 lj/sdk lj9_6 0.449100 3.731000 +pair_coeff 7 8 lj/sdk lj9_6 0.437200 4.033000 +pair_coeff 7 9 lj/sdk lj9_6 0.365000 3.987000 +pair_coeff 7 10 lj/sdk lj9_6 0.380000 3.840000 +pair_coeff 7 11 lj/sdk lj9_6 0.403000 3.981500 +pair_coeff 7 12 lj/sdk lj9_6 0.410000 3.890000 ### +pair_coeff 7 13 lj/sdk lj9_6 0.419100 3.731000 ### +pair_coeff 8 8 lj/sdk lj9_6 0.469000 4.585000 +pair_coeff 8 9 lj/sdk lj9_6 0.444000 4.545500 +pair_coeff 8 10 lj/sdk lj9_6 0.382500 4.403000 +pair_coeff 8 11 lj/sdk lj9_6 0.452700 4.165000 +pair_coeff 8 12 lj/sdk lj9_6 0.410000 4.340000 +pair_coeff 8 13 lj/sdk lj9_6 0.437200 4.033000 +pair_coeff 9 9 lj/sdk lj9_6 0.420000 4.506000 # +pair_coeff 9 10 lj/sdk lj9_6 0.362000 4.363500 # +pair_coeff 9 11 lj/sdk lj9_6 0.428400 4.378000 +pair_coeff 9 12 lj/sdk lj9_6 0.377000 4.274000 # +pair_coeff 9 13 lj/sdk lj9_6 0.365000 3.987000 # +pair_coeff 10 10 lj/sdk lj9_6 0.312000 4.221000 # +pair_coeff 10 11 lj/sdk lj9_6 0.369200 4.235500 +pair_coeff 10 12 lj/sdk lj9_6 0.370000 4.274000 # +pair_coeff 10 13 lj/sdk lj9_6 0.380000 3.840000 # +pair_coeff 11 11 lj/sdk lj9_6 0.437000 4.250000 +pair_coeff 11 12 lj/sdk lj9_6 0.420000 4.250000 +pair_coeff 11 13 lj/sdk lj9_6 0.443000 3.981500 +pair_coeff 12 12 lj/sdk lj9_6 0.405000 4.250000 # +pair_coeff 12 13 lj/sdk lj9_6 0.440000 3.890000 # +pair_coeff 13 13 lj/sdk lj9_6 0.449100 3.731000 # +pair_coeff 1 14 table tab.rc20.0_gold_water GOLD-WATER +pair_coeff 2 14 table tab.rc20.0_gold_CM GOLD-CM +pair_coeff 3 14 table tab.rc20.0_gold_CM GOLD-CM +pair_coeff 4 14 table tab.rc20.0_gold_CT2 GOLD-CT2 +pair_coeff 5 14 table tab.rc20.0_gold_EO GOLD-EO +pair_coeff 6 14 table tab.rc20.0_gold_EO GOLD-EO +pair_coeff 7 14 table tab.rc20.0_gold_OA GOLD-OA +pair_coeff 8 14 table tab.rc20.0_gold_OA GOLD-OA +pair_coeff 9 14 table tab.rc20.0_gold_CM GOLD-CM +pair_coeff 10 14 table tab.rc20.0_gold_CT2 GOLD-CT2 +pair_coeff 11 14 table tab.rc20.0_gold_EO GOLD-EO +pair_coeff 12 14 table tab.rc20.0_gold_EO GOLD-EO +pair_coeff 13 14 table tab.rc20.0_gold_OA GOLD-OA +pair_coeff 14 14 table tab.rc20.0_gold_OA GOLD-OA + + +bond_coeff 1 6.955000 3.710000 +bond_coeff 2 6.160000 3.640000 +bond_coeff 3 6.160000 3.650000 +bond_coeff 4 9.000000 3.130000 +bond_coeff 5 63.000000 2.160000 +bond_coeff 6 5.500000 3.330000 +bond_coeff 7 5.400000 3.340000 +bond_coeff 8 4.900000 3.280000 +bond_coeff 9 15.000000 2.790000 +bond_coeff 10 7.100000 3.610000 +bond_coeff 11 7.100000 3.560000 +bond_coeff 12 10.000000 3.070000 +bond_coeff 13 6.955000 3.050000 +bond_coeff 14 7.500000 3.010000 +bond_coeff 15 14.000000 2.530000 + +angle_coeff 1 1.190000 173.000000 lj9_6 0.420 4.5060 +angle_coeff 2 1.190000 175.000000 lj9_6 0.444 4.5455 +angle_coeff 3 1.600000 172.000000 lj9_6 0.362 4.3635 +angle_coeff 4 1.700000 173.000000 lj9_6 0.312 4.2210 +angle_coeff 5 1.093000 175.500000 lj9_6 0.469 4.5850 +angle_coeff 6 3.500000 135.000000 lj9_6 0.420 4.2500 +angle_coeff 7 3.400000 132.000000 lj9_6 0.405 4.2500 +angle_coeff 8 3.500000 134.000000 lj9_6 0.437 4.2500 +angle_coeff 9 6.600000 131.000000 lj9_6 0.4491 3.7310 +angle_coeff 10 3.000000 131.000000 lj9_6 0.440 3.8900 +angle_coeff 11 2.200000 145.000000 lj9_6 0.410 4.3400 +angle_coeff 12 1.300000 178.000000 lj9_6 0.370 4.1400 +angle_coeff 13 1.800000 171.000000 lj9_6 0.420 4.5060 +angle_coeff 14 2.600000 165.000000 lj9_6 0.312 4.2210 +angle_coeff 15 4.000000 146.000000 lj9_6 0.3800 3.8400 +angle_coeff 16 1.500000 172.000000 lj9_6 0.377 4.274 +angle_coeff 17 3.200000 146.000000 lj9_6 0.377 4.274 +angle_coeff 18 1.800000 166.000000 lj9_6 0.380 3.840 + diff --git a/C12E2_work/pteros_example.cxx b/C12E2_work/pteros_example.cxx new file mode 100644 index 00000000..5b5711f6 --- /dev/null +++ b/C12E2_work/pteros_example.cxx @@ -0,0 +1,23 @@ +#include +#include +#include +#include + +pteros::System load_system(std::string input_data) { + /* + sample system loading function + */ + pteros::System sample_system; + sample_system(input_data); + return sample_system; +} + + + +int main() { + pteros::System sys1; + sys1.load("/home/sang/Desktop/CG_NP_ANALYSIS/CG/NP/HYDROPHOBIC/1_NP_sim/1_NP_sim_NPT.gro"); + pteros::Selection sel0 = sys1("all"); + std::cout << sel0 << std::endl; + return 0; +} diff --git a/C12E2_work/valgrind-out.txt b/C12E2_work/valgrind-out.txt new file mode 100644 index 00000000..2b2b953b --- /dev/null +++ b/C12E2_work/valgrind-out.txt @@ -0,0 +1,20 @@ +==38755== Memcheck, a memory error detector +==38755== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al. +==38755== Using Valgrind-3.18.1-42b08ed5bd-20211015 and LibVEX; rerun with -h for copyright info +==38755== Command: ./new +==38755== Parent PID: 9232 +==38755== +--38755-- +--38755-- Valgrind options: +--38755-- --leak-check=full +--38755-- --show-leak-kinds=all +--38755-- --track-origins=yes +--38755-- --verbose +--38755-- --log-file=valgrind-out.txt +--38755-- Contents of /proc/version: +--38755-- Linux version 5.15.0-48-generic (buildd@lcy02-amd64-080) (gcc (Ubuntu 11.2.0-19ubuntu1) 11.2.0, GNU ld (GNU Binutils for Ubuntu) 2.38) #54-Ubuntu SMP Fri Aug 26 13:26:29 UTC 2022 +--38755-- +--38755-- Arch and hwcaps: AMD64, LittleEndian, amd64-cx16-rdtscp-sse3-ssse3-avx-f16c +--38755-- Page sizes: currently 4096, max supported 4096 +--38755-- Valgrind library directory: /usr/libexec/valgrind +--38755-- Reading syms from /home/sang/Desktop/git/Research/C12E2_work/new