NAMD
GlobalMasterTMD.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 "GlobalMasterTMD.h"
14 #include "PDB.h"
15 #include "PDBData.h"
16 #include "fitrms.h"
17 //#define DEBUGM
18 #define MIN_DEBUG_LEVEL 1
19 #include "Debug.h"
20 //using namespace __gnu_cxx;
21 using namespace std;
22 
23 class Matrix4TMD {
24  BigReal mat[16];
25 public:
26  Matrix4TMD(void) { identity(); }
27  Matrix4TMD(const BigReal *m) { memcpy(mat, m, 16*sizeof(BigReal)); }
28  void multpoint(BigReal point[3]) const {
29  BigReal tmp[3];
30  BigReal itmp3 = 1.0f / (point[0]*mat[3] + point[1]*mat[7] +
31  point[2]*mat[11] + mat[15]);
32  tmp[0] = itmp3*point[0];
33  tmp[1] = itmp3*point[1];
34  tmp[2] = itmp3*point[2];
35  point[0]=tmp[0]*mat[0] + tmp[1]*mat[4] + tmp[2]*mat[ 8] + itmp3*mat[12];
36  point[1]=tmp[0]*mat[1] + tmp[1]*mat[5] + tmp[2]*mat[ 9] + itmp3*mat[13];
37  point[2]=tmp[0]*mat[2] + tmp[1]*mat[6] + tmp[2]*mat[10] + itmp3*mat[14];
38  }
39 
40  void identity() {
41  memset(mat, 0, 16*sizeof(BigReal));
42  mat[0]=1.0f;
43  mat[5]=1.0f;
44  mat[10]=1.0f;
45  mat[15]=1.0f;
46  }
47  void transpose() {
48  BigReal tmp[16];
49  int i,j;
50  for(i=0;i<4;i++) {
51  for(j=0;j<4;j++) {
52  tmp[4*i+j] = mat[i+4*j];
53  }
54  }
55  for(i=0;i<16;i++) mat[i] = tmp[i];
56  }
58  void multmatrix(const Matrix4TMD &m) {
59  BigReal tmp[4];
60  for (int j=0; j<4; j++) {
61  tmp[0] = mat[j];
62  tmp[1] = mat[4+j];
63  tmp[2] = mat[8+j];
64  tmp[3] = mat[12+j];
65  for (int i=0; i<4; i++) {
66  mat[4*i+j] = m.mat[4*i]*tmp[0] + m.mat[4*i+1]*tmp[1] +
67  m.mat[4*i+2]*tmp[2] + m.mat[4*i+3]*tmp[3];
68  }
69  }
70  }
72  Matrix4TMD m;
73  m.mat[12] = x;
74  m.mat[13] = y;
75  m.mat[14] = z;
76  multmatrix(m);
77  }
78  void translate(BigReal d[3]) { translate(d[0], d[1], d[2]); }
79 };
80 
82  DebugM(3,"initialize called\n");
84  outputFreq = params->TMDOutputFreq;
85  K = params->TMDk;
86 /* if (params->TMDInitialRMSD < 0){
87  initialRMS = -1; // get from first coordinates
88  }
89  else */
90  initialRMS = params->TMDInitialRMSD;
91  finalRMS = params->TMDFinalRMSD;
92 
93  currentStep = params->firstTimestep;
94  firstStep = params->TMDFirstStep;
95  lastStep = params->TMDLastStep;
96  qDiffRMSD=params->TMDDiffRMSD;
97  altloc = 0;
98  target = 0;
99  target2 = 0;
100  weight = 0;
101  if (qDiffRMSD) parseAtoms(params->TMDFile2,Node::Object()->molecule->numAtoms, 1);
102  parseAtoms(params->TMDFile,Node::Object()->molecule->numAtoms, 0);
103 
104  //iterate through all domains to see if altloc is used
105  map <int, vector<int> >::iterator it;
106  for (it = dmap.begin(); it != dmap.end(); ++it){
107  int refcount = 0;
108  int biascount = 0;
109  for(int i = 0; i<it->second.size(); i++){
110  char aloc = altloc[it->second[i]];
111  if ( aloc & 1 ) ++biascount;
112  if ( aloc & 2 ) ++refcount;
113  }
114  altlocmap[it->first] = ( refcount ? 1 : 0 );
115  if ( ! refcount ) refcount = it->second.size();
116  iout << iINFO << "TMD domain " << it->first <<
117  " has " << it->second.size() << " atoms " <<
118  refcount << " fitted " << biascount << " biased\n" << endi;
119  }
120 
121  // k /= numTMDatoms;
122  iout << iINFO << numTMDatoms << " TMD ATOMS\n" << endi;
123  DebugM(1,"done with initialize\n");
124 }
125 
126 void GlobalMasterTMD::parseAtoms(const char *file, int numTotalAtoms, bool isTwo) {
127  DebugM(3,"parseAtoms called\n");
128  PDB tmdpdb(file);
129  numatoms = tmdpdb.num_atoms();
130  if (numatoms < 1)
131  NAMD_die("No atoms found in TMDFile\n");
132  if (numatoms < numTotalAtoms)
133  iout << iWARN << "The number of atoms in TMDFile is less than the total number of atoms in the structure.\n" << endi;
134  if (numatoms > numTotalAtoms)
135  NAMD_die("The number of atoms in TMDFile must not exceed the total number of atoms in the structure!");
136  if ( modifyRequestedAtoms().size() )
137  NAMD_bug("GlobalMasterTMD::parseAtoms() modifyRequestedAtoms() not empty");
138 
139  numTMDatoms = 0;
140  if(isTwo){
141  target2 = new BigReal[3*numatoms];
142  atompos2 = new Vector[numatoms];
143  tmdpdb.get_all_positions(atompos2);
144 
145  }
146  else{
147  target = new BigReal[3*numatoms];
148  atompos = new Vector[numatoms];
149  tmdpdb.get_all_positions(atompos);
150  }
151  if ( ! altloc ) altloc = new char[numatoms];
152  int i;
153  for (i=0; i<numatoms; i++) {
154 #ifdef MEM_OPT_VERSION
155  PDBCoreData *atom = tmdpdb.atom(i);
156  char aloc = tmdpdb.alternatelocation(i);
157 #else
158  PDBAtom *atom = tmdpdb.atom(i); // get an atom from the file
159  char aloc = atom->alternatelocation()[0];
160 #endif
161  if ( aloc ) aloc -= '0';
162  if ( aloc ) aloc = 2; // 2 bit == reference
163  if ( atom->occupancy() ) aloc |= 1; // 1 bit == biased
164  altloc[i] = aloc;
165  if ( aloc ) {
166  if(isTwo){
167  target2[3*numTMDatoms ] = atompos2[i].x;
168  target2[3*numTMDatoms+1] = atompos2[i].y;
169  target2[3*numTMDatoms+2] = atompos2[i].z;
170  }
171  else{
172  target[3*numTMDatoms ] = atompos[i].x;
173  target[3*numTMDatoms+1] = atompos[i].y;
174  target[3*numTMDatoms+2] = atompos[i].z;
175  // aidmap[i] = numTMDatoms++;
176  numTMDatoms++;
177  // add the atom to the list
178  modifyRequestedAtoms().add(i);
179  if(!K){ kmap[i] = atom->occupancy();}
180  //check to see if domain is already in the map
181  map <int, vector<int> >::iterator it = dmap.find((int)atom->temperaturefactor());
182  if (it != dmap.end()){
183  it->second.push_back(i); //add atomid to vector in proper existing domain
184  }
185  else{
186  dmap[(int)atom->temperaturefactor()] = vector <int> (1,i); //create new domain with atomid
187  }
188  }
189  }
190  }
191 
192  DebugM(1,"done with parseAtoms\n");
193 }
194 
195 //recreates target array for domain selected
196 void GlobalMasterTMD::NewTarget(int domain)
197 {
198  map <int, vector<int> >::iterator it = dmap.find(domain);
199  //target_aid = new int[it->second.size()];
200  delete [] target;
201  target = new BigReal [3*it->second.size()];
202  for(int i = 0; i<it->second.size(); i++){
203  target[3*i ] = atompos[it->second[i]].x;
204  target[3*i+1] = atompos[it->second[i]].y;
205  target[3*i+2] = atompos[it->second[i]].z;
206  }
207  if(qDiffRMSD){
208  delete [] target2;
209  target2 = new BigReal [3*it->second.size()];
210  for(int i = 0; i<it->second.size(); i++){
211  target2[3*i ] = atompos2[it->second[i]].x;
212  target2[3*i+1] = atompos2[it->second[i]].y;
213  target2[3*i+2] = atompos2[it->second[i]].z;
214  }
215  }
216  delete [] weight; weight = 0;
217  if ( altlocmap.find(domain)->second ) {
218  weight = new BigReal [it->second.size()];
219  for(int i = 0; i<it->second.size(); i++){
220  weight[i] = ( (altloc[it->second[i]] & 2) ? 1.0 : 0.0 );
221  }
222  }
223 // target_aid[i] = it->second[i];
224 // aidmap[it->second[i]] = i;
225 
226 }
228  delete [] target;
229  // delete [] target_aid;
230 // delete [] aidmap;
231  delete [] atompos;
232  delete [] atompos2;
233  delete [] altloc;
234  delete [] target2;
235 }
236 void GlobalMasterTMD::calculate() {
237  // have to reset my forces every time.
238  modifyAppliedForces().resize(0);
239  modifyForcedAtoms().resize(0);
240 
241  // see if TMD should be active
242  if (currentStep < firstStep || currentStep >= lastStep) {
243  currentStep++;
244  return;
245  }
246 
247  map<int, Position> positions;
248  AtomIDList::const_iterator a_i = getAtomIdBegin();
249  AtomIDList::const_iterator a_e = getAtomIdEnd();
250  PositionList::const_iterator p_i = getAtomPositionBegin();
251 
252  //create mapping of positions now to avoid
253  //going through these iterators for each domain
254  for ( ; a_i != a_e; ++a_i, ++p_i ){
255  positions[*a_i] = *p_i;
256  }
257 
258  //iterate through all domains
259  map <int, vector<int> >::iterator it;
260  for (it = dmap.begin(); it != dmap.end(); ++it){
261  NewTarget(it->first); //set new target
262  // fetch the current coordinates
263  BigReal *curpos = new BigReal[3*(it->second.size())];
264  for (int i = 0; i < it->second.size(); i++){
265  //int ind = 3*aidmap[it->second[i]];
266  curpos[3*i ] = positions[it->second[i]].x;
267  curpos[3*i+1] = positions[it->second[i]].y;
268  curpos[3*i+2] = positions[it->second[i]].z;
269  }
270  BigReal *curpos2;
271 if(qDiffRMSD){
272  curpos2 = new BigReal[3*(it->second.size())];
273  for (int i = 0; i < it->second.size(); i++){
274  //int ind = 3*aidmap[it->second[i]];
275  curpos2[3*i ] = positions[it->second[i]].x;
276  curpos2[3*i+1] = positions[it->second[i]].y;
277  curpos2[3*i+2] = positions[it->second[i]].z;
278  }
279 }
280  // align target with current coordinates. Uses same weight for all
281  // atoms. Maybe instead use weight from occupancy?
282  BigReal ttt[16], pre[3], post[3];
283  BigReal curRMS = MatrixFitRMS(it->second.size(), target, curpos, weight, ttt);
284  // Compute targetRMS.
285  if (initialRMS < 0) {
286  initialRMS = curRMS;
287  }
288 
289  BigReal curRMS0 = curRMS;
290  BigReal curRMS1 = 1.;
291  BigReal ttt1[16];
292  if(qDiffRMSD){
293  curRMS1 = MatrixFitRMS(it->second.size(), target2, curpos2, weight, ttt1);
294  curRMS = curRMS0 - curRMS1 ;
295  }
296 
297 
298  BigReal frac = (BigReal(currentStep-firstStep)) /
299  (lastStep-firstStep);
300 
301  BigReal targetRMS = initialRMS * (1-frac) + frac * finalRMS;
302 
303 
304  BigReal maxforce2 = 0.;
305 //orig finalRMS < initialRMS...changed to <= when allowing initialRMS = 0
306 //qdiff part and the whole && section new finalRMS <=initialRMS
307  if (((finalRMS <= initialRMS && targetRMS <= curRMS) ||
308  (finalRMS >= initialRMS && targetRMS >= curRMS) ||
309  qDiffRMSD) && (curRMS0 > 0. && curRMS1 > 0) ) {
310 
311 
312  // compute transformation to align target structure with current structure
313  // Might be more stable to instead align current positions with target,
314  // although then we have to back-transform the forces.
315  int j;
316  for (j=0; j<3; j++) {
317  post[j] = ttt[4*j+3];
318  ttt[4*j+3] = 0;
319  pre[j] = ttt[12+j];
320  ttt[12+j] = 0;
321  }
322  Matrix4TMD result;
323  result.translate(pre);
324  result.multmatrix(Matrix4TMD(ttt));
325  result.translate(post);
326 
327  // compute forces on each atom
328  BigReal myrms = 0;
329  for (int i=0; i<it->second.size(); i++) {
330  BigReal k = 0.;
331  if(!K){
332  k = kmap[it->second[i]];
333  } else if ( (! weight) || (altloc[it->second[i]] & 1) ) {
334  k = K/it->second.size();
335  }
336  // BigReal prefac = k * (targetRMS / curRMS - 1);
337  BigReal prefac = k * (targetRMS - curRMS)/curRMS0;
338  result.multpoint(target+3*i);
339  BigReal dx = curpos[3*i ] - target[3*i ];
340  BigReal dy = curpos[3*i+1] - target[3*i+1];
341  BigReal dz = curpos[3*i+2] - target[3*i+2];
342 
343  BigReal fvec[3] = { dx, dy, dz };
344  Vector force(fvec[0]*prefac, fvec[1]*prefac, fvec[2]*prefac);
345  modifyForcedAtoms().add(it->second[i]);
346  modifyAppliedForces().add(force);
347  BigReal force2 = force.length2();
348  if ( force2 > maxforce2 ) maxforce2 = force2;
349  }
350 
351  if(qDiffRMSD){
352  int j;
353  for (j=0; j<3; j++) {
354  post[j] = ttt1[4*j+3];
355  ttt1[4*j+3] = 0;
356  pre[j] = ttt1[12+j];
357  ttt1[12+j] = 0;
358  }
359  Matrix4TMD result2;
360  result2.identity();
361  result2.translate(pre);
362  result2.multmatrix(Matrix4TMD(ttt1));
363  result2.translate(post);
364 
365  // compute forces on each atom
366  BigReal myrms = 0;
367  for (int i=0; i<it->second.size(); i++) {
368  BigReal k = 0.;
369  if(!K){
370  k = kmap[it->second[i]];
371  } else if ( (! weight) || (altloc[it->second[i]] & 1) ) {
372  k = K/it->second.size();
373  }
374  // BigReal prefac = k * (targetRMS / curRMS - 1);
375  // BigReal prefac = k * (targetRMS - curRMS)/curRMS0;
376  BigReal prefac1 = - k * (targetRMS - curRMS)/curRMS1; // included with a negative term in the potential
377 
378  result2.multpoint(target2+3*i);
379  BigReal dx = curpos2[3*i ] - target2[3*i ];
380  BigReal dy = curpos2[3*i+1] - target2[3*i+1];
381  BigReal dz = curpos2[3*i+2] - target2[3*i+2];
382  BigReal fvec[3] = { dx, dy, dz };
383  Vector force(fvec[0]*prefac1, fvec[1]*prefac1, fvec[2]*prefac1);
384  modifyForcedAtoms().add(it->second[i]);
385  modifyAppliedForces().add(force);
386  BigReal force2 = force.length2();
387  if ( force2 > maxforce2 ) maxforce2 = force2;
388  }
389  }
390  }
391  //delete [] target_aid;
392  delete [] curpos;
393  if(qDiffRMSD){delete [] curpos2;}
394 // write output if needed
395  if (currentStep % outputFreq == 0) {
396  iout << "TMD " << currentStep << " Domain: "<< it->first << " " << targetRMS << ' ' << curRMS << '\n' << endi; //*it
397  // iout << "TMD " << currentStep << " " << targetRMS << ' ' << curRMS << ' ' << sqrt(maxforce2) << '\n' << endi;
398  }
399  }
400 /* // write output if needed
401  if (currentStep % outputFreq == 0) {
402  iout << "TMD " << currentStep << " " << targetRMS << ' ' << curRMS << '\n' << endi;
403  // iout << "TMD " << currentStep << " " << targetRMS << ' ' << curRMS << ' ' << sqrt(maxforce2) << '\n' << endi;
404  }*/
405  currentStep++;
406  // delete [] curpos;
407 }
408 
static Node * Object()
Definition: Node.h:86
void multmatrix(const Matrix4TMD &m)
premultiply the matrix by the given matrix, this-&gt;other * this
BigReal TMDInitialRMSD
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Definition: PDB.h:35
const Elem * const_iterator
Definition: ResizeArray.h:38
BigReal temperaturefactor(void)
Definition: PDBData.C:450
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
static int numatoms
Definition: ScriptTcl.C:64
#define DebugM(x, y)
Definition: Debug.h:59
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
char TMDFile2[128]
void transpose()
#define iout
Definition: InfoStream.h:87
Matrix4TMD(const BigReal *m)
void identity()
void NAMD_bug(const char *err_msg)
Definition: common.C:123
void translate(BigReal x, BigReal y, BigReal z)
gridSize z
const char * alternatelocation(void)
Definition: PDBData.C:392
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
BigReal MatrixFitRMS(int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
Definition: fitrms.C:91
BigReal TMDFinalRMSD
gridSize y
BigReal occupancy(void)
Definition: PDBData.C:444
infostream & endi(infostream &s)
Definition: InfoStream.C:38
gridSize x
char TMDFile[128]
Molecule * molecule
Definition: Node.h:176
void translate(BigReal d[3])
Matrix4TMD(void)
void multpoint(BigReal point[3]) const
double BigReal
Definition: common.h:112