GlobalMasterSymmetry.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "Node.h"
00009 #include "Molecule.h"
00010 #include "NamdTypes.h"
00011 
00012 #include "SimParameters.h"
00013 #include "GlobalMasterSymmetry.h"
00014 #include "PDB.h"
00015 #include "PDBData.h"
00016 #include "fitrms.h"
00017 #include <fstream>
00018 #include <sstream>
00019 #include <algorithm>
00020 #include <iterator>
00021 #include <stdlib.h>
00022 //#define DEBUGM
00023 #define MIN_DEBUG_LEVEL 1
00024 #include "Debug.h"
00025 //using namespace __gnu_cxx;
00026 using namespace std;
00027 //Read in and parse matrix file
00028 void GlobalMasterSymmetry::parseMatrix(int id, char fileName []){
00029   int count = 1;
00030   string line;
00031   ifstream matrixFile (fileName);
00032   if (! matrixFile.is_open()){
00033     NAMD_err((std::string("Error opening symmetry matrix file ") + fileName).c_str());
00034   } else {
00035     while (!matrixFile.eof()){
00036       vector <string> tmp;
00037       for (int i = 0; i < 4; i++){
00038       getline(matrixFile, line);
00039       istringstream iss(line);
00040       copy(istream_iterator<string>(iss),
00041             istream_iterator<string>(),
00042             back_inserter<vector <string>  >(tmp));
00043       }
00044       getline(matrixFile, line);
00045       Matrix4Symmetry tmpmatrix;
00046       if (tmp.size() < 16){NAMD_die("Error reading matrix file.  Please check layout of the matrix file(s).");}
00047       for(int j = 0; j < 16; j++){
00048         tmpmatrix.mat[j] = atof(tmp[j].c_str());
00049       }
00050       tmpmatrix.transpose();
00051       matrices.push_back(tmpmatrix);
00052       count++;
00053     }
00054     matrixFile.close();
00055   } 
00056 }
00057 GlobalMasterSymmetry::GlobalMasterSymmetry() {
00058   DebugM(3,"initialize called\n");
00059   SimParameters *params = Node::Object()->simParameters;
00060   currentStep = params->firstTimestep;
00061   firstStep = params->symmetryFirstStep;
00062   lastStep = params->symmetryLastStep;
00063   firstFullStep = params->symmetryFirstFullStep;
00064   lastFullStep = params->symmetryLastFullStep;
00065   K = params->symmetryk;
00066 
00067   StringList *klist = Node::Object()->configList->find("symmetrykfile");
00068   if (!K){
00069     //if (!params->symmetrykfile){NAMD_die("A pdb file containing per-atom force constants must be specified if symmetryk is not in configuration file!");}
00070     if (!klist){NAMD_die("A pdb file containing per-atom force constants must be specified if symmetryk is not in configuration file!");}
00071     //symmetrykfile = params->symmetrykfile;
00072   }
00073   scaleForces = params->symmetryScaleForces;
00074   if (scaleForces && lastStep == -1){
00075     NAMD_die("symmetryLastStep must be specified if symmetryScaleForces is enabled!");
00076   }
00077   StringList *matrixList = Node::Object()->configList->find("symmetryMatrixFile");
00078   StringList *symmetryList = Node::Object()->configList->find("symmetryFile");
00079   int symfileindex = 0;
00080   if (!K) {symmetrykfile = klist->data;}
00081   for (; symmetryList; symmetryList = symmetryList->next) {  
00082     parseAtoms(symmetryList->data, Node::Object()->molecule->numAtoms, symfileindex);
00083     if(!K){
00084       klist = klist->next;
00085      if (klist){ symmetrykfile = klist->data;}
00086     }
00087   }
00088 
00089   map<int, vector<int> >::iterator it = simmap.begin();
00090   if (!matrixList) {initialTransform();}
00091   for (; matrixList; matrixList = matrixList->next, ++it){
00092       parseMatrix(it->first, matrixList->data);
00093   }
00094 
00095   DebugM(1,"done with initialize\n");
00096 }
00097 
00098 //Aligns monomers based on transformation matrices
00099 //found in matrix file
00100 /*
00101 void GlobalMasterSymmetry::alignMonomers(){
00102   //this is assuming the matrices are written
00103   //in order of monomer id designation (0, 1, 2,..etc)
00104   map<int, vector<int> >::iterator simit = simmap.begin();
00105   for (; simit!=simmap.end(); ++simit){
00106     for(int x = 0; x < simit->second.size(); x++){
00107     map<int, vector<int> >::iterator mit = dmap.find(simit->second[x]);
00108     for (int i = 0; i < mit->second.size(); i++){
00109       map<int, BigReal *>::iterator it = posmap.find(mit->second[i]);
00110       matrices[(mit->first)-1].multpoint(it->second);
00111     }
00112   }
00113   }
00114 }
00115 */
00116 bool GlobalMasterSymmetry::gluInvertMatrix(const BigReal m[16], BigReal invOut[16])
00117 {
00118   BigReal inv[16], det;
00119   int i;
00120 
00121   inv[0] =   m[5]*m[10]*m[15] - m[5]*m[11]*m[14] - m[9]*m[6]*m[15]
00122   + m[9]*m[7]*m[14] + m[13]*m[6]*m[11] - m[13]*m[7]*m[10];
00123   inv[4] =  -m[4]*m[10]*m[15] + m[4]*m[11]*m[14] + m[8]*m[6]*m[15]
00124   - m[8]*m[7]*m[14] - m[12]*m[6]*m[11] + m[12]*m[7]*m[10];
00125   inv[8] =   m[4]*m[9]*m[15] - m[4]*m[11]*m[13] - m[8]*m[5]*m[15]
00126   + m[8]*m[7]*m[13] + m[12]*m[5]*m[11] - m[12]*m[7]*m[9];
00127   inv[12] = -m[4]*m[9]*m[14] + m[4]*m[10]*m[13] + m[8]*m[5]*m[14]
00128   - m[8]*m[6]*m[13] - m[12]*m[5]*m[10] + m[12]*m[6]*m[9];
00129   inv[1] =  -m[1]*m[10]*m[15] + m[1]*m[11]*m[14] + m[9]*m[2]*m[15]
00130   - m[9]*m[3]*m[14] - m[13]*m[2]*m[11] + m[13]*m[3]*m[10];
00131   inv[5] =   m[0]*m[10]*m[15] - m[0]*m[11]*m[14] - m[8]*m[2]*m[15]
00132   + m[8]*m[3]*m[14] + m[12]*m[2]*m[11] - m[12]*m[3]*m[10];
00133   inv[9] =  -m[0]*m[9]*m[15] + m[0]*m[11]*m[13] + m[8]*m[1]*m[15]
00134   - m[8]*m[3]*m[13] - m[12]*m[1]*m[11] + m[12]*m[3]*m[9];
00135   inv[13] =  m[0]*m[9]*m[14] - m[0]*m[10]*m[13] - m[8]*m[1]*m[14]
00136   + m[8]*m[2]*m[13] + m[12]*m[1]*m[10] - m[12]*m[2]*m[9];
00137   inv[2] =   m[1]*m[6]*m[15] - m[1]*m[7]*m[14] - m[5]*m[2]*m[15]
00138   + m[5]*m[3]*m[14] + m[13]*m[2]*m[7] - m[13]*m[3]*m[6];
00139   inv[6] =  -m[0]*m[6]*m[15] + m[0]*m[7]*m[14] + m[4]*m[2]*m[15]
00140   - m[4]*m[3]*m[14] - m[12]*m[2]*m[7] + m[12]*m[3]*m[6];
00141   inv[10] =  m[0]*m[5]*m[15] - m[0]*m[7]*m[13] - m[4]*m[1]*m[15]
00142   + m[4]*m[3]*m[13] + m[12]*m[1]*m[7] - m[12]*m[3]*m[5];
00143   inv[14] = -m[0]*m[5]*m[14] + m[0]*m[6]*m[13] + m[4]*m[1]*m[14]
00144   - m[4]*m[2]*m[13] - m[12]*m[1]*m[6] + m[12]*m[2]*m[5];
00145   inv[3] =  -m[1]*m[6]*m[11] + m[1]*m[7]*m[10] + m[5]*m[2]*m[11]
00146   - m[5]*m[3]*m[10] - m[9]*m[2]*m[7] + m[9]*m[3]*m[6];
00147   inv[7] =   m[0]*m[6]*m[11] - m[0]*m[7]*m[10] - m[4]*m[2]*m[11]
00148   + m[4]*m[3]*m[10] + m[8]*m[2]*m[7] - m[8]*m[3]*m[6];
00149   inv[11] = -m[0]*m[5]*m[11] + m[0]*m[7]*m[9] + m[4]*m[1]*m[11]
00150   - m[4]*m[3]*m[9] - m[8]*m[1]*m[7] + m[8]*m[3]*m[5];
00151   inv[15] =  m[0]*m[5]*m[10] - m[0]*m[6]*m[9] - m[4]*m[1]*m[10]
00152   + m[4]*m[2]*m[9] + m[8]*m[1]*m[6] - m[8]*m[2]*m[5];
00153 
00154   det = m[0]*inv[0] + m[1]*inv[4] + m[2]*inv[8] + m[3]*inv[12];
00155   if (det == 0)
00156             return false;
00157 
00158             det = 1.0 / det;
00159 
00160             for (i = 0; i < 16; i++)
00161                       invOut[i] = inv[i] * det;
00162 
00163                       return true;
00164 }
00165 
00166 void GlobalMasterSymmetry::initialTransform(){
00167   map<int, vector<int> >::iterator simit = simmap.begin();
00168 
00169   for (; simit!=simmap.end(); ++simit){
00170   map<int, vector<int> >::iterator fit = dmap.find(simit->second[0]);
00171   BigReal * first = new BigReal [3*fit->second.size()];
00172   for(int i = 0; i < fit->second.size(); i++){
00173     map<int, BigReal *>::iterator sit = posmap.find(fit->second[i]);
00174     first[3*i] = sit->second[0];
00175     first[3*i+1] = sit->second[1];
00176     first[3*i+2] = sit->second[2];  
00177   }
00178 
00179   map<int, vector<int> >::iterator it = dmap.begin();
00180   for(; it != dmap.end(); ++it){
00181   if (std::find(simit->second.begin(), simit->second.end(), it->first)!=simit->second.end()){
00182   BigReal * arr = new BigReal [3*it->second.size()];
00183     for (int i = 0; i<it->second.size(); i++){
00184       map<int, BigReal *>::iterator sit = posmap.find(it->second[i]);
00185       arr[3*i] = sit->second[0];
00186       arr[3*i+1] = sit->second[1];
00187       arr[3*i+2] = sit->second[2];  
00188     }  
00189       BigReal ttt[16], pre[3], post[3];
00190       BigReal curRMS = MatrixFitRMS(it->second.size(), arr, first, NULL, ttt);
00191       int j;
00192       for (j=0; j<3; j++) {
00193         post[j] = ttt[4*j+3];
00194         ttt[4*j+3] = 0;
00195         pre[j] = ttt[12+j];
00196         ttt[12+j] = 0;
00197       }
00198       Matrix4Symmetry result;
00199       result.translate(pre);
00200       result.multmatrix(Matrix4Symmetry(ttt));
00201       result.translate(post);
00202       matrices.push_back(result);
00203 
00204       delete [] arr;
00205   }
00206   }
00207   delete [] first;
00208   }
00209 
00210   for (int i = 0; i < matrices.size(); i++) {
00211     matrices[i].transpose();
00212     //iout << "Matrix: " << i << " " << matrices[i].mat[0] << " " << matrices[i].mat[1] << " " << matrices[i].mat[2] << " "<<matrices[i].mat[3] << "\n";
00213     //iout << "Matrix: " << i << " " << matrices[i].mat[4] << " " << matrices[i].mat[5] << " " << matrices[i].mat[6] << " "<<matrices[i].mat[7] << "\n";  
00214     //iout << "Matrix: " << i << " " << matrices[i].mat[8] << " " << matrices[i].mat[9] << " " << matrices[i].mat[10] << " "<<matrices[i].mat[11] << "\n";
00215     //iout << "Matrix: " << i << " " << matrices[i].mat[12] << " " << matrices[i].mat[13] << " " << matrices[i].mat[14] << " "<<matrices[i].mat[15] << "\n";
00216     //iout <<"\n"<<endi;
00217     matrices[i].transpose();
00218   }
00219 }
00220 
00221 void GlobalMasterSymmetry::parseAtoms(const char *file, int numTotalAtoms, int symfileindex) {
00222   DebugM(3,"parseAtoms called\n");
00223   PDB tmdpdb(file);
00224   int numatoms = tmdpdb.num_atoms();
00225   if (numatoms < 1) 
00226     NAMD_die("No atoms found in symmetryFile\n");
00227   if (numatoms < numTotalAtoms)
00228     iout << iWARN << "The number of atoms in symmetryFile is less than the total number of atoms in the structure.\n" << endi;
00229   if (numatoms > numTotalAtoms)
00230     NAMD_die("The number of atoms in symmetryFile must not exceed the total number of atoms in the structure!");
00231  // if ( modifyRequestedAtoms().size() )
00232  //   NAMD_bug("GlobalMasterSymmetry::parseAtoms() modifyRequestedAtoms() not empty");
00233 
00234   Vector * atompos = new Vector[numatoms];  
00235   tmdpdb.get_all_positions(atompos);
00236   if (K){symmetrykfile = file;} 
00237   PDB kpdb(symmetrykfile);
00238 
00239   int i;
00240   for (i=0; i<numatoms; i++) {
00241 #ifdef MEM_OPT_VERSION
00242     PDBCoreData *atom = tmdpdb.atom(i);
00243 #else
00244     PDBAtom *atom = tmdpdb.atom(i); // get an atom from the file
00245 #endif
00246 
00247     if (atom->occupancy() && atom->temperaturefactor()) { // if occupancy and beta are not 0, then add it!
00248       // add the atom to the list
00249       modifyRequestedAtoms().add(i);
00250       if(!K){
00251         #ifdef MEM_OPT_VERSION
00252           PDBCoreData *atomk = kpdb.atom(i);
00253         #else
00254           PDBAtom *atomk = kpdb.atom(i); // get an atom from the file
00255         #endif 
00256         //kmap[i] = atomk->occupancy();
00257         kdmap[atom->temperaturefactor()][i] = atomk->occupancy();
00258       }
00259       BigReal *arr = new BigReal [3];
00260       arr[0] = atompos[i].x;
00261       arr[1] = atompos[i].y;
00262       arr[2] = atompos[i].z;
00263       posmap[i] = arr;
00264 
00265       bmap[atom->temperaturefactor()] = atom->occupancy();
00266         //check to see if monomer id is already in the map
00267       map <int, vector<int> >::iterator it = dmap.find((int)atom->temperaturefactor());
00268       map <int, vector<int> >::iterator simit = simmap.find((int)atom->occupancy());
00269       if (it != dmap.end()){
00270         it->second.push_back(i); //add atomid to vector in proper existing monomer id
00271       }
00272       else{
00273          dmap[(int)atom->temperaturefactor()] = vector <int> (1,i); //create new monomer id with atomid
00274       }
00275       if (simit != simmap.end()){
00276         if (std::find(simit->second.begin(), simit->second.end(), atom->temperaturefactor()) == simit->second.end()){
00277         simit->second.push_back(atom->temperaturefactor()); 
00278         }
00279       }
00280       else {
00281         simmap[(int)atom->occupancy()] = vector <int> (1, atom->temperaturefactor());
00282       }
00283     }
00284   }
00285 
00286    map <int, vector<int> >::iterator simit = simmap.begin();
00287    for (; simit != simmap.end(); ++simit){
00288     map <int, vector<int> >::iterator sit = dmap.find(simit->second[0]);
00289     int numatoms = sit->second.size();
00290     for (int i = 0; i<simit->second.size(); i++){
00291       map <int, vector<int> >::iterator fit = dmap.find(simit->second[i]);
00292       if (fit->second.size() != numatoms){
00293         NAMD_die("Every monomer must contain the same number of atoms!");
00294       }
00295     }
00296    }  
00297   delete [] atompos;
00298 }
00299 
00300 void GlobalMasterSymmetry::determineAverage() {
00301 
00302    map <int, vector<BigReal *> >::iterator delit = averagePos.begin();
00303    for (; delit != averagePos.end(); ++delit){
00304      for (int i = 0; i < delit->second.size(); i++){
00305        delete [] delit->second[i];
00306      }
00307      delit->second.erase(delit->second.begin(), delit->second.end());
00308    }
00309    //std::map <int, BigReal * > posmap;
00310    map <int, BigReal *>::iterator posit;
00311    map <int, vector<int> >::iterator simit = simmap.begin();
00312    for (; simit != simmap.end(); ++simit){     
00313     
00314 
00315     map <int, BigReal *>::iterator pit = posmap.begin();
00316     for (; pit != posmap.end(); ++pit){delete [] pit->second;}
00317     posmap.clear();
00318 
00319     map <int, vector<int> >::iterator dit = dmap.begin();
00320     for (; dit!=dmap.end(); ++dit){
00321       for (int i = 0; i < dit->second.size(); i++){
00322         if (std::find(simit->second.begin(), simit->second.end(), dit->first)!=simit->second.end()){
00323 
00324           BigReal* arr = new BigReal[3];
00325           arr[0] = positions[dit->second[i]].x;
00326           arr[1] = positions[dit->second[i]].y;
00327           arr[2] = positions[dit->second[i]].z;
00328 
00329           posmap[dit->second[i]] = arr;
00330         }
00331       } 
00332     }
00333     averagePos[simit->first] = vector <BigReal *> (); 
00334     map <int, vector<int> >::iterator it = dmap.begin();
00335     map <int, vector<int> >::iterator sit = dmap.find(simit->second[0]);
00336     int numatoms = sit->second.size();
00337     for (int i = 0; i < numatoms; i++){
00338        BigReal *arr = new BigReal [3];
00339       arr[0] = 0;
00340       arr[1] = 0;
00341       arr[2] = 0;
00342       for (; it!=dmap.end(); ++it){
00343         if (std::find(simit->second.begin(), simit->second.end(), it->first)!=simit->second.end()){
00344           posit = posmap.find(it->second[i]);
00345           matrices[(it->first)-1].multpoint(posit->second);
00346           arr[0] += posit->second[0];
00347           arr[1] += posit->second[1];
00348           arr[2] += posit->second[2];
00349         }
00350       }
00351       it = dmap.begin();
00352       BigReal *avg = new BigReal[3];
00353       avg[0] = arr[0]/(simit->second.size());
00354       avg[1] = arr[1]/(simit->second.size());
00355       avg[2] = arr[2]/(simit->second.size());
00356       averagePos[simit->first].push_back(avg);
00357       delete [] arr;
00358     }
00359     
00360    }
00361 
00362 }
00363 
00364 void GlobalMasterSymmetry::backTransform(){
00365   map <int, BigReal *>::iterator bit = backavg.begin();
00366   for (; bit != backavg.end(); ++bit){delete [] bit->second;}
00367   backavg.clear();
00368 
00369   map <int, vector<int> >::iterator it = dmap.begin();
00370   for (; it!=dmap.end(); ++it){
00371     map<int, int >::iterator bmit = bmap.find(it->first);
00372     int bm = bmit->second;
00373     map<int, vector <BigReal *> >::iterator avit = averagePos.find(bmit->second);
00374     int numatoms = it->second.size();
00375     BigReal *avg = new BigReal [3*numatoms];
00376     for (int i = 0; i < numatoms; i++){
00377       avg[3*i] = avit->second[i][0];
00378       avg[3*i+1] = avit->second[i][1];
00379       avg[3*i+2] = avit->second[i][2];
00380     }
00381     BigReal inverse[16];
00382     matrices[it->first-1].transpose();
00383     gluInvertMatrix(matrices[it->first-1].mat, inverse);
00384     Matrix4Symmetry inv(inverse);
00385     inv.transpose();
00386     matrices[it->first-1].transpose();
00387  
00388     for (int k = 0; k < numatoms; k++){
00389       inv.multpoint(avg+3*k);
00390     }
00391     backavg[it->first] = avg;
00392    }
00393 
00394 }
00395 GlobalMasterSymmetry::~GlobalMasterSymmetry() { 
00396   map <int, BigReal *>::iterator pit = posmap.begin();
00397   for (; pit != posmap.end(); ++pit){delete pit->second;}
00398 }
00399 void GlobalMasterSymmetry::calculate() {
00400   // have to reset my forces every time.  
00401   modifyAppliedForces().resize(0);
00402   modifyForcedAtoms().resize(0);
00403 
00404   // see if symmetry restraints should be active
00405   if (currentStep < firstStep || (currentStep >= lastStep && lastStep != -1)) {
00406     currentStep++;
00407     return;
00408   }
00409 
00410   //map<int, Position> positions;
00411   AtomIDList::const_iterator a_i = getAtomIdBegin();
00412   AtomIDList::const_iterator a_e = getAtomIdEnd();
00413   PositionList::const_iterator p_i = getAtomPositionBegin();
00414 
00415   //create mapping of positions now to avoid
00416   //going through these iterators for each domain
00417   for ( ; a_i != a_e; ++a_i, ++p_i ){
00418     positions[*a_i] = *p_i;
00419   }
00420 
00421   map <int, vector<int> >::iterator it;
00422   //map <int, BigReal *>::iterator pit = posmap.begin();
00423   //for (; pit != posmap.end(); ++pit){delete [] pit->second;}
00424 
00425  // posmap.clear();
00426   //for (it = dmap.begin(); it != dmap.end(); ++it){
00427     // fetch the current coordinates
00428     
00429     //for (int i = 0; i < it->second.size(); i++){
00430 
00431       //BigReal* arr = new BigReal[3];
00432       //arr[0] = positions[it->second[i]].x;
00433       //arr[1] = positions[it->second[i]].y;
00434       //arr[2] = positions[it->second[i]].z;
00435 
00436      // posmap[it->second[i]] = arr;
00437    // } 
00438 //}
00439 
00440 //  alignMonomers();
00441   determineAverage();
00442   backTransform();
00443 
00444   //iterate through all domains
00445   for (it = dmap.begin(); it != dmap.end(); ++it){
00446 
00447   // fetch the current coordinates
00448   BigReal *curpos = new BigReal[3*(it->second.size())];
00449   for (int i = 0; i < it->second.size(); i++){
00450     curpos[3*i  ] = positions[it->second[i]].x;
00451     curpos[3*i+1] = positions[it->second[i]].y;
00452     curpos[3*i+2] = positions[it->second[i]].z;  
00453   }
00454 
00455 
00456   BigReal *tmpavg = backavg[it->first];
00457 
00458   for (int i=0; i<it->second.size(); i++) {
00459     BigReal k; 
00460     if(!K){
00461      //k = kmap[it->second[i]];  
00462      k = kdmap[it->first][it->second[i]];
00463     }
00464     else{
00465      k = K/it->second.size();  
00466     }
00467     BigReal maxforce2 = 0.;
00468     BigReal frac = -k;
00469     if (scaleForces){
00470 
00471     if (currentStep < firstFullStep){
00472       BigReal linear_evolve = (BigReal(currentStep-firstStep)) /
00473                             (firstFullStep-firstStep);
00474       frac = frac*(linear_evolve);
00475     }
00476     if (currentStep > lastFullStep)
00477     {
00478       BigReal linear_evolve = (BigReal(currentStep-lastFullStep)) /
00479                             (lastStep-lastFullStep);
00480       frac = frac*(1-linear_evolve);
00481     }
00482 
00483     }
00484     BigReal dx = curpos[3*i  ] - tmpavg[3*i];
00485 //    iout << "DX: " << dx << "\n" << endi;
00486     BigReal dy = curpos[3*i+1] - tmpavg[3*i+1];
00487     BigReal dz = curpos[3*i+2] - tmpavg[3*i+2];
00488     BigReal fvec[3] = { dx, dy, dz };
00489     Vector force(fvec[0]*frac, fvec[1]*frac, fvec[2]*frac);
00490     modifyForcedAtoms().add(it->second[i]);
00491     modifyAppliedForces().add(force);
00492     BigReal force2 = force.length2();
00493     if ( force2 > maxforce2 ) maxforce2 = force2;
00494   } 
00495   delete [] curpos;
00496  }
00497    
00498   currentStep++;
00499 }

Generated on Wed Nov 22 01:17:15 2017 for NAMD by  doxygen 1.4.7