GlobalMasterTMD.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 "GlobalMasterTMD.h"
00014 #include "PDB.h"
00015 #include "PDBData.h"
00016 #include "fitrms.h"
00017 //#define DEBUGM
00018 #define MIN_DEBUG_LEVEL 1
00019 #include "Debug.h"
00020 //using namespace __gnu_cxx;
00021 using namespace std;
00022 
00023 class Matrix4TMD {
00024   BigReal mat[16];                               
00025 public:
00026   Matrix4TMD(void) { identity(); }
00027   Matrix4TMD(const BigReal *m)  { memcpy(mat, m, 16*sizeof(BigReal)); }
00028   void multpoint(BigReal point[3]) const {
00029     BigReal tmp[3];
00030     BigReal itmp3 = 1.0f / (point[0]*mat[3] + point[1]*mat[7] +
00031                             point[2]*mat[11] + mat[15]);
00032     tmp[0] = itmp3*point[0];
00033     tmp[1] = itmp3*point[1];
00034     tmp[2] = itmp3*point[2];
00035     point[0]=tmp[0]*mat[0] + tmp[1]*mat[4] + tmp[2]*mat[ 8] + itmp3*mat[12];
00036     point[1]=tmp[0]*mat[1] + tmp[1]*mat[5] + tmp[2]*mat[ 9] + itmp3*mat[13];
00037     point[2]=tmp[0]*mat[2] + tmp[1]*mat[6] + tmp[2]*mat[10] + itmp3*mat[14];
00038   }
00039 
00040   void identity() {
00041     memset(mat, 0, 16*sizeof(BigReal));
00042     mat[0]=1.0f;
00043     mat[5]=1.0f;
00044     mat[10]=1.0f;
00045     mat[15]=1.0f;
00046   }
00047   void transpose() {
00048     BigReal tmp[16];
00049     int i,j;
00050     for(i=0;i<4;i++) {
00051       for(j=0;j<4;j++) {
00052         tmp[4*i+j] = mat[i+4*j];
00053       }
00054     }
00055     for(i=0;i<16;i++) mat[i] = tmp[i];
00056   }
00058   void multmatrix(const Matrix4TMD &m) {
00059     BigReal tmp[4];
00060     for (int j=0; j<4; j++) {
00061       tmp[0] = mat[j];
00062       tmp[1] = mat[4+j];
00063       tmp[2] = mat[8+j]; 
00064       tmp[3] = mat[12+j];
00065       for (int i=0; i<4; i++) {
00066         mat[4*i+j] = m.mat[4*i]*tmp[0] + m.mat[4*i+1]*tmp[1] +
00067           m.mat[4*i+2]*tmp[2] + m.mat[4*i+3]*tmp[3]; 
00068       }
00069     } 
00070   }
00071   void translate(BigReal x, BigReal y, BigReal z) {
00072     Matrix4TMD m;               
00073     m.mat[12] = x;
00074     m.mat[13] = y;
00075     m.mat[14] = z;
00076     multmatrix(m);
00077   }
00078   void translate(BigReal d[3]) { translate(d[0], d[1], d[2]); }
00079 };
00080 
00081 GlobalMasterTMD::GlobalMasterTMD() {
00082   DebugM(3,"initialize called\n");
00083   SimParameters *params = Node::Object()->simParameters;
00084   outputFreq = params->TMDOutputFreq;
00085   K = params->TMDk;
00086 /*  if (params->TMDInitialRMSD < 0){
00087     initialRMS = -1; // get from first coordinates
00088   }
00089   else */
00090     initialRMS = params->TMDInitialRMSD;
00091   finalRMS = params->TMDFinalRMSD;
00092   
00093   currentStep = params->firstTimestep;
00094   firstStep = params->TMDFirstStep;
00095   lastStep = params->TMDLastStep;
00096   qDiffRMSD=params->TMDDiffRMSD;
00097   altloc = 0;
00098   target = 0;
00099   target2 = 0;
00100   weight = 0;
00101   if (qDiffRMSD) parseAtoms(params->TMDFile2,Node::Object()->molecule->numAtoms, 1);
00102   parseAtoms(params->TMDFile,Node::Object()->molecule->numAtoms, 0);
00103 
00104   //iterate through all domains to see if altloc is used
00105   map <int, vector<int> >::iterator it;
00106   for (it = dmap.begin(); it != dmap.end(); ++it){  
00107     int refcount = 0;
00108     int biascount = 0;
00109     for(int i = 0; i<it->second.size(); i++){
00110       char aloc = altloc[it->second[i]];
00111       if ( aloc & 1 ) ++biascount;
00112       if ( aloc & 2 ) ++refcount;
00113     }
00114     altlocmap[it->first] = ( refcount ? 1 : 0 );
00115     if ( ! refcount ) refcount = it->second.size();
00116     iout << iINFO << "TMD domain " << it->first <<
00117       " has " << it->second.size() << " atoms " <<
00118       refcount << " fitted " << biascount << " biased\n" << endi;
00119   }
00120 
00121  // k /= numTMDatoms;
00122   iout << iINFO << numTMDatoms << " TMD ATOMS\n" << endi;
00123   DebugM(1,"done with initialize\n");
00124 }
00125 
00126 void GlobalMasterTMD::parseAtoms(const char *file, int numTotalAtoms, bool isTwo) {
00127   DebugM(3,"parseAtoms called\n");
00128   PDB tmdpdb(file);
00129   numatoms = tmdpdb.num_atoms();
00130   if (numatoms < 1) 
00131     NAMD_die("No atoms found in TMDFile\n");
00132   if (numatoms < numTotalAtoms)
00133     iout << iWARN << "The number of atoms in TMDFile is less than the total number of atoms in the structure.\n" << endi;
00134   if (numatoms > numTotalAtoms)
00135     NAMD_die("The number of atoms in TMDFile must not exceed the total number of atoms in the structure!");
00136   if ( modifyRequestedAtoms().size() )
00137     NAMD_bug("GlobalMasterTMD::parseAtoms() modifyRequestedAtoms() not empty");
00138 
00139   numTMDatoms = 0;
00140   if(isTwo){
00141     target2 = new BigReal[3*numatoms];
00142     atompos2 = new Vector[numatoms];
00143     tmdpdb.get_all_positions(atompos2);
00144 
00145   }
00146   else{
00147     target = new BigReal[3*numatoms];
00148     atompos = new Vector[numatoms];
00149     tmdpdb.get_all_positions(atompos);
00150   }
00151   if ( ! altloc ) altloc = new char[numatoms];
00152   int i;
00153   for (i=0; i<numatoms; i++) {
00154 #ifdef MEM_OPT_VERSION
00155     PDBCoreData *atom = tmdpdb.atom(i);
00156     char aloc = tmdpdb.alternatelocation(i);
00157 #else
00158     PDBAtom *atom = tmdpdb.atom(i); // get an atom from the file
00159     char aloc = atom->alternatelocation()[0];
00160 #endif
00161     if ( aloc ) aloc -= '0';
00162     if ( aloc ) aloc = 2;  // 2 bit == reference
00163     if ( atom->occupancy() ) aloc |= 1;  // 1 bit == biased
00164     altloc[i] = aloc;
00165     if ( aloc ) {
00166       if(isTwo){
00167         target2[3*numTMDatoms  ] = atompos2[i].x;
00168         target2[3*numTMDatoms+1] = atompos2[i].y;
00169         target2[3*numTMDatoms+2] = atompos2[i].z;
00170       }
00171       else{
00172         target[3*numTMDatoms  ] = atompos[i].x;
00173         target[3*numTMDatoms+1] = atompos[i].y;
00174         target[3*numTMDatoms+2] = atompos[i].z;
00175  //       aidmap[i] = numTMDatoms++;
00176         numTMDatoms++;
00177         // add the atom to the list
00178         modifyRequestedAtoms().add(i);
00179         if(!K){ kmap[i] = atom->occupancy();}
00180         //check to see if domain is already in the map
00181         map <int, vector<int> >::iterator it = dmap.find((int)atom->temperaturefactor());
00182         if (it != dmap.end()){
00183           it->second.push_back(i); //add atomid to vector in proper existing domain
00184         }
00185         else{
00186            dmap[(int)atom->temperaturefactor()] = vector <int> (1,i); //create new domain with atomid
00187         }
00188       }
00189     }
00190   }
00191   
00192   DebugM(1,"done with parseAtoms\n");
00193 }
00194 
00195 //recreates target array for domain selected
00196 void GlobalMasterTMD::NewTarget(int domain)
00197 {
00198   map <int, vector<int> >::iterator it = dmap.find(domain);
00199   //target_aid = new int[it->second.size()];
00200   delete [] target;
00201   target = new BigReal [3*it->second.size()];
00202   for(int i = 0; i<it->second.size(); i++){
00203    target[3*i  ] = atompos[it->second[i]].x;
00204    target[3*i+1] = atompos[it->second[i]].y;
00205    target[3*i+2] = atompos[it->second[i]].z; 
00206   }
00207   if(qDiffRMSD){
00208    delete [] target2;
00209    target2 = new BigReal [3*it->second.size()];
00210    for(int i = 0; i<it->second.size(); i++){
00211      target2[3*i  ] = atompos2[it->second[i]].x;
00212      target2[3*i+1] = atompos2[it->second[i]].y;
00213      target2[3*i+2] = atompos2[it->second[i]].z; 
00214    }   
00215   }
00216   delete [] weight;  weight = 0;
00217   if ( altlocmap.find(domain)->second ) {
00218     weight = new BigReal [it->second.size()];
00219     for(int i = 0; i<it->second.size(); i++){
00220       weight[i] = ( (altloc[it->second[i]] & 2) ? 1.0 : 0.0 );
00221     }
00222   }
00223 //   target_aid[i] = it->second[i];
00224 //   aidmap[it->second[i]] = i; 
00225 
00226 }
00227 GlobalMasterTMD::~GlobalMasterTMD() { 
00228   delete [] target;
00229  // delete [] target_aid;
00230 //  delete [] aidmap;
00231   delete [] atompos;
00232   delete [] atompos2;
00233   delete [] altloc;
00234   delete [] target2;
00235 }
00236 void GlobalMasterTMD::calculate() {
00237   // have to reset my forces every time.  
00238   modifyAppliedForces().resize(0);
00239   modifyForcedAtoms().resize(0);
00240   
00241   // see if TMD should be active
00242   if (currentStep < firstStep || currentStep >= lastStep) {
00243     currentStep++;
00244     return;
00245   }
00246   
00247   map<int, Position> positions;
00248   AtomIDList::const_iterator a_i = getAtomIdBegin();
00249   AtomIDList::const_iterator a_e = getAtomIdEnd();
00250   PositionList::const_iterator p_i = getAtomPositionBegin();
00251 
00252   //create mapping of positions now to avoid
00253   //going through these iterators for each domain
00254   for ( ; a_i != a_e; ++a_i, ++p_i ){
00255     positions[*a_i] = *p_i;
00256   }
00257 
00258   //iterate through all domains
00259   map <int, vector<int> >::iterator it;
00260   for (it = dmap.begin(); it != dmap.end(); ++it){  
00261   NewTarget(it->first);  //set new target
00262   // fetch the current coordinates
00263   BigReal *curpos = new BigReal[3*(it->second.size())];
00264   for (int i = 0; i < it->second.size(); i++){
00265     //int ind = 3*aidmap[it->second[i]];
00266     curpos[3*i  ] = positions[it->second[i]].x;
00267     curpos[3*i+1] = positions[it->second[i]].y;
00268     curpos[3*i+2] = positions[it->second[i]].z;
00269   }
00270   BigReal *curpos2;
00271 if(qDiffRMSD){
00272   curpos2 = new BigReal[3*(it->second.size())];
00273   for (int i = 0; i < it->second.size(); i++){
00274     //int ind = 3*aidmap[it->second[i]];
00275     curpos2[3*i  ] = positions[it->second[i]].x;
00276     curpos2[3*i+1] = positions[it->second[i]].y;
00277     curpos2[3*i+2] = positions[it->second[i]].z;
00278   }
00279 }
00280   // align target with current coordinates.   Uses same weight for all
00281   // atoms.  Maybe instead use weight from occupancy?
00282   BigReal ttt[16], pre[3], post[3];
00283   BigReal curRMS = MatrixFitRMS(it->second.size(), target, curpos, weight, ttt);
00284   // Compute targetRMS.
00285   if (initialRMS < 0) {
00286     initialRMS = curRMS;
00287   }
00288 
00289   BigReal curRMS0 = curRMS;
00290   BigReal curRMS1 =  1.; 
00291   BigReal ttt1[16]; 
00292   if(qDiffRMSD){
00293     curRMS1 = MatrixFitRMS(it->second.size(), target2, curpos2, weight, ttt1);
00294     curRMS = curRMS0 - curRMS1 ;
00295   }
00296 
00297 
00298   BigReal frac = (BigReal(currentStep-firstStep)) /
00299                             (lastStep-firstStep);
00300   
00301   BigReal targetRMS = initialRMS * (1-frac) + frac * finalRMS;
00302 
00303   
00304   BigReal maxforce2 = 0.;
00305 //orig finalRMS < initialRMS...changed to <= when allowing initialRMS = 0
00306 //qdiff part and the whole && section new finalRMS <=initialRMS
00307   if (((finalRMS <= initialRMS && targetRMS <= curRMS) ||
00308       (finalRMS >= initialRMS && targetRMS >= curRMS) ||
00309       qDiffRMSD) && (curRMS0 > 0. && curRMS1 > 0) ) {
00310 
00311 
00312     // compute transformation to align target structure with current structure
00313     // Might be more stable to instead align current positions with target,
00314     // although then we have to back-transform the forces.
00315     int j;
00316     for (j=0; j<3; j++) {
00317       post[j] = ttt[4*j+3];
00318       ttt[4*j+3] = 0;
00319       pre[j] = ttt[12+j];
00320       ttt[12+j] = 0;
00321     }
00322     Matrix4TMD result;
00323     result.translate(pre);
00324     result.multmatrix(Matrix4TMD(ttt));
00325     result.translate(post);
00326   
00327     // compute forces on each atom
00328     BigReal myrms = 0;
00329       for (int i=0; i<it->second.size(); i++) {
00330       BigReal k = 0.; 
00331       if(!K){
00332         k = kmap[it->second[i]];  
00333       } else if ( (! weight) || (altloc[it->second[i]] & 1) ) {
00334         k = K/it->second.size();  
00335       }
00336    //   BigReal prefac = k * (targetRMS / curRMS - 1); 
00337       BigReal prefac = k * (targetRMS - curRMS)/curRMS0;
00338       result.multpoint(target+3*i);
00339       BigReal dx = curpos[3*i  ] - target[3*i  ];
00340       BigReal dy = curpos[3*i+1] - target[3*i+1];
00341       BigReal dz = curpos[3*i+2] - target[3*i+2];
00342       
00343       BigReal fvec[3] = { dx, dy, dz };
00344       Vector force(fvec[0]*prefac, fvec[1]*prefac, fvec[2]*prefac);
00345       modifyForcedAtoms().add(it->second[i]);
00346       modifyAppliedForces().add(force);
00347       BigReal force2 = force.length2();
00348       if ( force2 > maxforce2 ) maxforce2 = force2;
00349     }
00350 
00351     if(qDiffRMSD){
00352        int j;
00353       for (j=0; j<3; j++) {
00354         post[j] = ttt1[4*j+3];
00355         ttt1[4*j+3] = 0;
00356         pre[j] = ttt1[12+j];
00357         ttt1[12+j] = 0;
00358       }
00359       Matrix4TMD result2;
00360       result2.identity();
00361       result2.translate(pre);
00362       result2.multmatrix(Matrix4TMD(ttt1));
00363       result2.translate(post);
00364     
00365       // compute forces on each atom
00366       BigReal myrms = 0;
00367         for (int i=0; i<it->second.size(); i++) {
00368         BigReal k = 0.;
00369         if(!K){
00370           k = kmap[it->second[i]];  
00371         } else if ( (! weight) || (altloc[it->second[i]] & 1) ) {
00372           k = K/it->second.size();  
00373         }
00374      //   BigReal prefac = k * (targetRMS / curRMS - 1); 
00375   //      BigReal prefac = k * (targetRMS - curRMS)/curRMS0;
00376           BigReal prefac1 = - k * (targetRMS  - curRMS)/curRMS1; // included with a negative term in the potential
00377 
00378         result2.multpoint(target2+3*i);
00379         BigReal dx = curpos2[3*i  ] - target2[3*i  ];
00380         BigReal dy = curpos2[3*i+1] - target2[3*i+1];
00381         BigReal dz = curpos2[3*i+2] - target2[3*i+2];
00382         BigReal fvec[3] = { dx, dy, dz };
00383         Vector force(fvec[0]*prefac1, fvec[1]*prefac1, fvec[2]*prefac1);
00384         modifyForcedAtoms().add(it->second[i]);
00385         modifyAppliedForces().add(force);
00386         BigReal force2 = force.length2();
00387         if ( force2 > maxforce2 ) maxforce2 = force2;
00388       }   
00389     }
00390   }
00391   //delete [] target_aid;
00392   delete [] curpos;
00393   if(qDiffRMSD){delete [] curpos2;}
00394 // write output if needed
00395   if (currentStep % outputFreq == 0) {
00396     iout << "TMD  " << currentStep << " Domain: "<< it->first << " " << targetRMS << ' ' << curRMS << '\n' << endi; //*it
00397     // iout << "TMD  " << currentStep << " " << targetRMS << ' ' << curRMS << ' ' << sqrt(maxforce2) << '\n' << endi;
00398   }
00399   }
00400 /*  // write output if needed
00401   if (currentStep % outputFreq == 0) {
00402     iout << "TMD  " << currentStep << " " << targetRMS << ' ' << curRMS << '\n' << endi;
00403     // iout << "TMD  " << currentStep << " " << targetRMS << ' ' << curRMS << ' ' << sqrt(maxforce2) << '\n' << endi;
00404   }*/
00405   currentStep++;
00406  // delete [] curpos;
00407 }
00408 

Generated on Tue Nov 21 01:17:13 2017 for NAMD by  doxygen 1.4.7