NAMD
GlobalMasterSymmetry.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "Node.h"
9 #include "Molecule.h"
10 #include "NamdTypes.h"
11 
12 #include "SimParameters.h"
13 #include "GlobalMasterSymmetry.h"
14 #include "PDB.h"
15 #include "PDBData.h"
16 #include "fitrms.h"
17 #include <fstream>
18 #include <sstream>
19 #include <algorithm>
20 #include <iterator>
21 #include <stdlib.h>
22 //#define DEBUGM
23 #define MIN_DEBUG_LEVEL 1
24 #include "Debug.h"
25 //using namespace __gnu_cxx;
26 using namespace std;
27 //Read in and parse matrix file
28 void GlobalMasterSymmetry::parseMatrix(int id, char fileName []){
29  int count = 1;
30  string line;
31  ifstream matrixFile (fileName);
32  if (! matrixFile.is_open()){
33  NAMD_err((std::string("Error opening symmetry matrix file ") + fileName).c_str());
34  } else {
35  while (!matrixFile.eof()){
36  vector <string> tmp;
37  for (int i = 0; i < 4; i++){
38  getline(matrixFile, line);
39  istringstream iss(line);
40  copy(istream_iterator<string>(iss),
41  istream_iterator<string>(),
42  back_inserter<vector <string> >(tmp));
43  }
44  getline(matrixFile, line);
45  Matrix4Symmetry tmpmatrix;
46  if (tmp.size() < 16){NAMD_die("Error reading matrix file. Please check layout of the matrix file(s).");}
47  for(int j = 0; j < 16; j++){
48  tmpmatrix.mat[j] = atof(tmp[j].c_str());
49  }
50  tmpmatrix.transpose();
51  matrices.push_back(tmpmatrix);
52  count++;
53  }
54  matrixFile.close();
55  }
56 }
58  DebugM(3,"initialize called\n");
60  currentStep = params->firstTimestep;
61  firstStep = params->symmetryFirstStep;
62  lastStep = params->symmetryLastStep;
63  firstFullStep = params->symmetryFirstFullStep;
64  lastFullStep = params->symmetryLastFullStep;
65  K = params->symmetryk;
66 
67  StringList *klist = Node::Object()->configList->find("symmetrykfile");
68  if (!K){
69  //if (!params->symmetrykfile){NAMD_die("A pdb file containing per-atom force constants must be specified if symmetryk is not in configuration file!");}
70  if (!klist){NAMD_die("A pdb file containing per-atom force constants must be specified if symmetryk is not in configuration file!");}
71  //symmetrykfile = params->symmetrykfile;
72  }
73  scaleForces = params->symmetryScaleForces;
74  if (scaleForces && lastStep == -1){
75  NAMD_die("symmetryLastStep must be specified if symmetryScaleForces is enabled!");
76  }
77  StringList *matrixList = Node::Object()->configList->find("symmetryMatrixFile");
78  StringList *symmetryList = Node::Object()->configList->find("symmetryFile");
79  int symfileindex = 0;
80  if (!K) {symmetrykfile = klist->data;}
81  for (; symmetryList; symmetryList = symmetryList->next) {
82  parseAtoms(symmetryList->data, Node::Object()->molecule->numAtoms, symfileindex);
83  if(!K){
84  klist = klist->next;
85  if (klist){ symmetrykfile = klist->data;}
86  }
87  }
88 
89  map<int, vector<int> >::iterator it = simmap.begin();
90  if (!matrixList) {initialTransform();}
91  for (; matrixList; matrixList = matrixList->next, ++it){
92  parseMatrix(it->first, matrixList->data);
93  }
94 
95  DebugM(1,"done with initialize\n");
96 }
97 
98 //Aligns monomers based on transformation matrices
99 //found in matrix file
100 /*
101 void GlobalMasterSymmetry::alignMonomers(){
102  //this is assuming the matrices are written
103  //in order of monomer id designation (0, 1, 2,..etc)
104  map<int, vector<int> >::iterator simit = simmap.begin();
105  for (; simit!=simmap.end(); ++simit){
106  for(int x = 0; x < simit->second.size(); x++){
107  map<int, vector<int> >::iterator mit = dmap.find(simit->second[x]);
108  for (int i = 0; i < mit->second.size(); i++){
109  map<int, BigReal *>::iterator it = posmap.find(mit->second[i]);
110  matrices[(mit->first)-1].multpoint(it->second);
111  }
112  }
113  }
114 }
115 */
116 bool GlobalMasterSymmetry::gluInvertMatrix(const BigReal m[16], BigReal invOut[16])
117 {
118  BigReal inv[16], det;
119  int i;
120 
121  inv[0] = m[5]*m[10]*m[15] - m[5]*m[11]*m[14] - m[9]*m[6]*m[15]
122  + m[9]*m[7]*m[14] + m[13]*m[6]*m[11] - m[13]*m[7]*m[10];
123  inv[4] = -m[4]*m[10]*m[15] + m[4]*m[11]*m[14] + m[8]*m[6]*m[15]
124  - m[8]*m[7]*m[14] - m[12]*m[6]*m[11] + m[12]*m[7]*m[10];
125  inv[8] = m[4]*m[9]*m[15] - m[4]*m[11]*m[13] - m[8]*m[5]*m[15]
126  + m[8]*m[7]*m[13] + m[12]*m[5]*m[11] - m[12]*m[7]*m[9];
127  inv[12] = -m[4]*m[9]*m[14] + m[4]*m[10]*m[13] + m[8]*m[5]*m[14]
128  - m[8]*m[6]*m[13] - m[12]*m[5]*m[10] + m[12]*m[6]*m[9];
129  inv[1] = -m[1]*m[10]*m[15] + m[1]*m[11]*m[14] + m[9]*m[2]*m[15]
130  - m[9]*m[3]*m[14] - m[13]*m[2]*m[11] + m[13]*m[3]*m[10];
131  inv[5] = m[0]*m[10]*m[15] - m[0]*m[11]*m[14] - m[8]*m[2]*m[15]
132  + m[8]*m[3]*m[14] + m[12]*m[2]*m[11] - m[12]*m[3]*m[10];
133  inv[9] = -m[0]*m[9]*m[15] + m[0]*m[11]*m[13] + m[8]*m[1]*m[15]
134  - m[8]*m[3]*m[13] - m[12]*m[1]*m[11] + m[12]*m[3]*m[9];
135  inv[13] = m[0]*m[9]*m[14] - m[0]*m[10]*m[13] - m[8]*m[1]*m[14]
136  + m[8]*m[2]*m[13] + m[12]*m[1]*m[10] - m[12]*m[2]*m[9];
137  inv[2] = m[1]*m[6]*m[15] - m[1]*m[7]*m[14] - m[5]*m[2]*m[15]
138  + m[5]*m[3]*m[14] + m[13]*m[2]*m[7] - m[13]*m[3]*m[6];
139  inv[6] = -m[0]*m[6]*m[15] + m[0]*m[7]*m[14] + m[4]*m[2]*m[15]
140  - m[4]*m[3]*m[14] - m[12]*m[2]*m[7] + m[12]*m[3]*m[6];
141  inv[10] = m[0]*m[5]*m[15] - m[0]*m[7]*m[13] - m[4]*m[1]*m[15]
142  + m[4]*m[3]*m[13] + m[12]*m[1]*m[7] - m[12]*m[3]*m[5];
143  inv[14] = -m[0]*m[5]*m[14] + m[0]*m[6]*m[13] + m[4]*m[1]*m[14]
144  - m[4]*m[2]*m[13] - m[12]*m[1]*m[6] + m[12]*m[2]*m[5];
145  inv[3] = -m[1]*m[6]*m[11] + m[1]*m[7]*m[10] + m[5]*m[2]*m[11]
146  - m[5]*m[3]*m[10] - m[9]*m[2]*m[7] + m[9]*m[3]*m[6];
147  inv[7] = m[0]*m[6]*m[11] - m[0]*m[7]*m[10] - m[4]*m[2]*m[11]
148  + m[4]*m[3]*m[10] + m[8]*m[2]*m[7] - m[8]*m[3]*m[6];
149  inv[11] = -m[0]*m[5]*m[11] + m[0]*m[7]*m[9] + m[4]*m[1]*m[11]
150  - m[4]*m[3]*m[9] - m[8]*m[1]*m[7] + m[8]*m[3]*m[5];
151  inv[15] = m[0]*m[5]*m[10] - m[0]*m[6]*m[9] - m[4]*m[1]*m[10]
152  + m[4]*m[2]*m[9] + m[8]*m[1]*m[6] - m[8]*m[2]*m[5];
153 
154  det = m[0]*inv[0] + m[1]*inv[4] + m[2]*inv[8] + m[3]*inv[12];
155  if (det == 0)
156  return false;
157 
158  det = 1.0 / det;
159 
160  for (i = 0; i < 16; i++)
161  invOut[i] = inv[i] * det;
162 
163  return true;
164 }
165 
166 void GlobalMasterSymmetry::initialTransform(){
167  map<int, vector<int> >::iterator simit = simmap.begin();
168 
169  for (; simit!=simmap.end(); ++simit){
170  map<int, vector<int> >::iterator fit = dmap.find(simit->second[0]);
171  BigReal * first = new BigReal [3*fit->second.size()];
172  for(int i = 0; i < fit->second.size(); i++){
173  map<int, BigReal *>::iterator sit = posmap.find(fit->second[i]);
174  first[3*i] = sit->second[0];
175  first[3*i+1] = sit->second[1];
176  first[3*i+2] = sit->second[2];
177  }
178 
179  map<int, vector<int> >::iterator it = dmap.begin();
180  for(; it != dmap.end(); ++it){
181  if (std::find(simit->second.begin(), simit->second.end(), it->first)!=simit->second.end()){
182  BigReal * arr = new BigReal [3*it->second.size()];
183  for (int i = 0; i<it->second.size(); i++){
184  map<int, BigReal *>::iterator sit = posmap.find(it->second[i]);
185  arr[3*i] = sit->second[0];
186  arr[3*i+1] = sit->second[1];
187  arr[3*i+2] = sit->second[2];
188  }
189  BigReal ttt[16], pre[3], post[3];
190  BigReal curRMS = MatrixFitRMS(it->second.size(), arr, first, NULL, ttt);
191  int j;
192  for (j=0; j<3; j++) {
193  post[j] = ttt[4*j+3];
194  ttt[4*j+3] = 0;
195  pre[j] = ttt[12+j];
196  ttt[12+j] = 0;
197  }
198  Matrix4Symmetry result;
199  result.translate(pre);
200  result.multmatrix(Matrix4Symmetry(ttt));
201  result.translate(post);
202  matrices.push_back(result);
203 
204  delete [] arr;
205  }
206  }
207  delete [] first;
208  }
209 
210  for (int i = 0; i < matrices.size(); i++) {
211  matrices[i].transpose();
212  //iout << "Matrix: " << i << " " << matrices[i].mat[0] << " " << matrices[i].mat[1] << " " << matrices[i].mat[2] << " "<<matrices[i].mat[3] << "\n";
213  //iout << "Matrix: " << i << " " << matrices[i].mat[4] << " " << matrices[i].mat[5] << " " << matrices[i].mat[6] << " "<<matrices[i].mat[7] << "\n";
214  //iout << "Matrix: " << i << " " << matrices[i].mat[8] << " " << matrices[i].mat[9] << " " << matrices[i].mat[10] << " "<<matrices[i].mat[11] << "\n";
215  //iout << "Matrix: " << i << " " << matrices[i].mat[12] << " " << matrices[i].mat[13] << " " << matrices[i].mat[14] << " "<<matrices[i].mat[15] << "\n";
216  //iout <<"\n"<<endi;
217  matrices[i].transpose();
218  }
219 }
220 
221 void GlobalMasterSymmetry::parseAtoms(const char *file, int numTotalAtoms, int symfileindex) {
222  DebugM(3,"parseAtoms called\n");
223  PDB tmdpdb(file);
224  int numatoms = tmdpdb.num_atoms();
225  if (numatoms < 1)
226  NAMD_die("No atoms found in symmetryFile\n");
227  if (numatoms < numTotalAtoms)
228  iout << iWARN << "The number of atoms in symmetryFile is less than the total number of atoms in the structure.\n" << endi;
229  if (numatoms > numTotalAtoms)
230  NAMD_die("The number of atoms in symmetryFile must not exceed the total number of atoms in the structure!");
231  // if ( modifyRequestedAtoms().size() )
232  // NAMD_bug("GlobalMasterSymmetry::parseAtoms() modifyRequestedAtoms() not empty");
233 
234  Vector * atompos = new Vector[numatoms];
235  tmdpdb.get_all_positions(atompos);
236  if (K){symmetrykfile = file;}
237  PDB kpdb(symmetrykfile);
238 
239  int i;
240  for (i=0; i<numatoms; i++) {
241 #ifdef MEM_OPT_VERSION
242  PDBCoreData *atom = tmdpdb.atom(i);
243 #else
244  PDBAtom *atom = tmdpdb.atom(i); // get an atom from the file
245 #endif
246 
247  if (atom->occupancy() && atom->temperaturefactor()) { // if occupancy and beta are not 0, then add it!
248  // add the atom to the list
249  modifyRequestedAtoms().add(i);
250  if(!K){
251  #ifdef MEM_OPT_VERSION
252  PDBCoreData *atomk = kpdb.atom(i);
253  #else
254  PDBAtom *atomk = kpdb.atom(i); // get an atom from the file
255  #endif
256  //kmap[i] = atomk->occupancy();
257  kdmap[atom->temperaturefactor()][i] = atomk->occupancy();
258  }
259  BigReal *arr = new BigReal [3];
260  arr[0] = atompos[i].x;
261  arr[1] = atompos[i].y;
262  arr[2] = atompos[i].z;
263  posmap[i] = arr;
264 
265  bmap[atom->temperaturefactor()] = atom->occupancy();
266  //check to see if monomer id is already in the map
267  map <int, vector<int> >::iterator it = dmap.find((int)atom->temperaturefactor());
268  map <int, vector<int> >::iterator simit = simmap.find((int)atom->occupancy());
269  if (it != dmap.end()){
270  it->second.push_back(i); //add atomid to vector in proper existing monomer id
271  }
272  else{
273  dmap[(int)atom->temperaturefactor()] = vector <int> (1,i); //create new monomer id with atomid
274  }
275  if (simit != simmap.end()){
276  if (std::find(simit->second.begin(), simit->second.end(), atom->temperaturefactor()) == simit->second.end()){
277  simit->second.push_back(atom->temperaturefactor());
278  }
279  }
280  else {
281  simmap[(int)atom->occupancy()] = vector <int> (1, atom->temperaturefactor());
282  }
283  }
284  }
285 
286  map <int, vector<int> >::iterator simit = simmap.begin();
287  for (; simit != simmap.end(); ++simit){
288  map <int, vector<int> >::iterator sit = dmap.find(simit->second[0]);
289  int numatoms = sit->second.size();
290  for (int i = 0; i<simit->second.size(); i++){
291  map <int, vector<int> >::iterator fit = dmap.find(simit->second[i]);
292  if (fit->second.size() != numatoms){
293  NAMD_die("Every monomer must contain the same number of atoms!");
294  }
295  }
296  }
297  delete [] atompos;
298 }
299 
300 void GlobalMasterSymmetry::determineAverage() {
301 
302  map <int, vector<BigReal *> >::iterator delit = averagePos.begin();
303  for (; delit != averagePos.end(); ++delit){
304  for (int i = 0; i < delit->second.size(); i++){
305  delete [] delit->second[i];
306  }
307  delit->second.erase(delit->second.begin(), delit->second.end());
308  }
309  //std::map <int, BigReal * > posmap;
310  map <int, BigReal *>::iterator posit;
311  map <int, vector<int> >::iterator simit = simmap.begin();
312  for (; simit != simmap.end(); ++simit){
313 
314 
315  map <int, BigReal *>::iterator pit = posmap.begin();
316  for (; pit != posmap.end(); ++pit){delete [] pit->second;}
317  posmap.clear();
318 
319  map <int, vector<int> >::iterator dit = dmap.begin();
320  for (; dit!=dmap.end(); ++dit){
321  for (int i = 0; i < dit->second.size(); i++){
322  if (std::find(simit->second.begin(), simit->second.end(), dit->first)!=simit->second.end()){
323 
324  BigReal* arr = new BigReal[3];
325  arr[0] = positions[dit->second[i]].x;
326  arr[1] = positions[dit->second[i]].y;
327  arr[2] = positions[dit->second[i]].z;
328 
329  posmap[dit->second[i]] = arr;
330  }
331  }
332  }
333  averagePos[simit->first] = vector <BigReal *> ();
334  map <int, vector<int> >::iterator it = dmap.begin();
335  map <int, vector<int> >::iterator sit = dmap.find(simit->second[0]);
336  int numatoms = sit->second.size();
337  for (int i = 0; i < numatoms; i++){
338  BigReal *arr = new BigReal [3];
339  arr[0] = 0;
340  arr[1] = 0;
341  arr[2] = 0;
342  for (; it!=dmap.end(); ++it){
343  if (std::find(simit->second.begin(), simit->second.end(), it->first)!=simit->second.end()){
344  posit = posmap.find(it->second[i]);
345  matrices[(it->first)-1].multpoint(posit->second);
346  arr[0] += posit->second[0];
347  arr[1] += posit->second[1];
348  arr[2] += posit->second[2];
349  }
350  }
351  it = dmap.begin();
352  BigReal *avg = new BigReal[3];
353  avg[0] = arr[0]/(simit->second.size());
354  avg[1] = arr[1]/(simit->second.size());
355  avg[2] = arr[2]/(simit->second.size());
356  averagePos[simit->first].push_back(avg);
357  delete [] arr;
358  }
359 
360  }
361 
362 }
363 
364 void GlobalMasterSymmetry::backTransform(){
365  map <int, BigReal *>::iterator bit = backavg.begin();
366  for (; bit != backavg.end(); ++bit){delete [] bit->second;}
367  backavg.clear();
368 
369  map <int, vector<int> >::iterator it = dmap.begin();
370  for (; it!=dmap.end(); ++it){
371  map<int, int >::iterator bmit = bmap.find(it->first);
372  int bm = bmit->second;
373  map<int, vector <BigReal *> >::iterator avit = averagePos.find(bmit->second);
374  int numatoms = it->second.size();
375  BigReal *avg = new BigReal [3*numatoms];
376  for (int i = 0; i < numatoms; i++){
377  avg[3*i] = avit->second[i][0];
378  avg[3*i+1] = avit->second[i][1];
379  avg[3*i+2] = avit->second[i][2];
380  }
381  BigReal inverse[16];
382  matrices[it->first-1].transpose();
383  gluInvertMatrix(matrices[it->first-1].mat, inverse);
384  Matrix4Symmetry inv(inverse);
385  inv.transpose();
386  matrices[it->first-1].transpose();
387 
388  for (int k = 0; k < numatoms; k++){
389  inv.multpoint(avg+3*k);
390  }
391  backavg[it->first] = avg;
392  }
393 
394 }
396  map <int, BigReal *>::iterator pit = posmap.begin();
397  for (; pit != posmap.end(); ++pit){delete pit->second;}
398 }
399 void GlobalMasterSymmetry::calculate() {
400  // have to reset my forces every time.
401  modifyAppliedForces().resize(0);
402  modifyForcedAtoms().resize(0);
403 
404  // see if symmetry restraints should be active
405  if (currentStep < firstStep || (currentStep >= lastStep && lastStep != -1)) {
406  currentStep++;
407  return;
408  }
409 
410  //map<int, Position> positions;
411  AtomIDList::const_iterator a_i = getAtomIdBegin();
412  AtomIDList::const_iterator a_e = getAtomIdEnd();
413  PositionList::const_iterator p_i = getAtomPositionBegin();
414 
415  //create mapping of positions now to avoid
416  //going through these iterators for each domain
417  for ( ; a_i != a_e; ++a_i, ++p_i ){
418  positions[*a_i] = *p_i;
419  }
420 
421  map <int, vector<int> >::iterator it;
422  //map <int, BigReal *>::iterator pit = posmap.begin();
423  //for (; pit != posmap.end(); ++pit){delete [] pit->second;}
424 
425  // posmap.clear();
426  //for (it = dmap.begin(); it != dmap.end(); ++it){
427  // fetch the current coordinates
428 
429  //for (int i = 0; i < it->second.size(); i++){
430 
431  //BigReal* arr = new BigReal[3];
432  //arr[0] = positions[it->second[i]].x;
433  //arr[1] = positions[it->second[i]].y;
434  //arr[2] = positions[it->second[i]].z;
435 
436  // posmap[it->second[i]] = arr;
437  // }
438 //}
439 
440 // alignMonomers();
441  determineAverage();
442  backTransform();
443 
444  //iterate through all domains
445  for (it = dmap.begin(); it != dmap.end(); ++it){
446 
447  // fetch the current coordinates
448  BigReal *curpos = new BigReal[3*(it->second.size())];
449  for (int i = 0; i < it->second.size(); i++){
450  curpos[3*i ] = positions[it->second[i]].x;
451  curpos[3*i+1] = positions[it->second[i]].y;
452  curpos[3*i+2] = positions[it->second[i]].z;
453  }
454 
455 
456  BigReal *tmpavg = backavg[it->first];
457 
458  for (int i=0; i<it->second.size(); i++) {
459  BigReal k;
460  if(!K){
461  //k = kmap[it->second[i]];
462  k = kdmap[it->first][it->second[i]];
463  }
464  else{
465  k = K/it->second.size();
466  }
467  BigReal maxforce2 = 0.;
468  BigReal frac = -k;
469  if (scaleForces){
470 
471  if (currentStep < firstFullStep){
472  BigReal linear_evolve = (BigReal(currentStep-firstStep)) /
473  (firstFullStep-firstStep);
474  frac = frac*(linear_evolve);
475  }
476  if (currentStep > lastFullStep)
477  {
478  BigReal linear_evolve = (BigReal(currentStep-lastFullStep)) /
479  (lastStep-lastFullStep);
480  frac = frac*(1-linear_evolve);
481  }
482 
483  }
484  BigReal dx = curpos[3*i ] - tmpavg[3*i];
485 // iout << "DX: " << dx << "\n" << endi;
486  BigReal dy = curpos[3*i+1] - tmpavg[3*i+1];
487  BigReal dz = curpos[3*i+2] - tmpavg[3*i+2];
488  BigReal fvec[3] = { dx, dy, dz };
489  Vector force(fvec[0]*frac, fvec[1]*frac, fvec[2]*frac);
490  modifyForcedAtoms().add(it->second[i]);
491  modifyAppliedForces().add(force);
492  BigReal force2 = force.length2();
493  if ( force2 > maxforce2 ) maxforce2 = force2;
494  }
495  delete [] curpos;
496  }
497 
498  currentStep++;
499 }
static Node * Object()
Definition: Node.h:86
Definition: PDB.h:35
const Elem * const_iterator
Definition: ResizeArray.h:38
void NAMD_err(const char *err_msg)
Definition: common.C:102
void multmatrix(const Matrix4Symmetry &)
premultiply the matrix by the given matrix, this-&gt;other * this
BigReal temperaturefactor(void)
Definition: PDBData.C:450
void translate(BigReal x, BigReal y, BigReal z)
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
BigReal mat[16]
static int numatoms
Definition: ScriptTcl.C:64
int symmetryLastFullStep
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
#define iout
Definition: InfoStream.h:87
Bool symmetryScaleForces
BigReal x
Definition: Vector.h:66
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
ConfigList * configList
Definition: Node.h:179
BigReal MatrixFitRMS(int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
Definition: fitrms.C:91
StringList * next
Definition: ConfigList.h:49
char * data
Definition: ConfigList.h:48
BigReal y
Definition: Vector.h:66
BigReal symmetryk
StringList * find(const char *name) const
Definition: ConfigList.C:341
int symmetryFirstFullStep
BigReal occupancy(void)
Definition: PDBData.C:444
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:112