23 #define MIN_DEBUG_LEVEL 1
28 void GlobalMasterSymmetry::parseMatrix(
int id,
char fileName []){
31 ifstream matrixFile (fileName);
32 if (! matrixFile.is_open()){
33 NAMD_err((std::string(
"Error opening symmetry matrix file ") + fileName).c_str());
35 while (!matrixFile.eof()){
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));
44 getline(matrixFile, line);
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());
51 matrices.push_back(tmpmatrix);
58 DebugM(3,
"initialize called\n");
70 if (!klist){
NAMD_die(
"A pdb file containing per-atom force constants must be specified if symmetryk is not in configuration file!");}
74 if (scaleForces && lastStep == -1){
75 NAMD_die(
"symmetryLastStep must be specified if symmetryScaleForces is enabled!");
80 if (!K) {symmetrykfile = klist->
data;}
81 for (; symmetryList; symmetryList = symmetryList->
next) {
85 if (klist){ symmetrykfile = klist->
data;}
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);
95 DebugM(1,
"done with initialize\n");
116 bool GlobalMasterSymmetry::gluInvertMatrix(
const BigReal m[16],
BigReal invOut[16])
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];
154 det = m[0]*inv[0] + m[1]*inv[4] + m[2]*inv[8] + m[3]*inv[12];
160 for (i = 0; i < 16; i++)
161 invOut[i] = inv[i] * det;
166 void GlobalMasterSymmetry::initialTransform(){
167 map<int, vector<int> >::iterator simit = simmap.begin();
169 for (; simit!=simmap.end(); ++simit){
170 map<int, vector<int> >::iterator fit = dmap.find(simit->second[0]);
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];
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()){
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];
189 BigReal ttt[16], pre[3], post[3];
192 for (j=0; j<3; j++) {
193 post[j] = ttt[4*j+3];
202 matrices.push_back(result);
210 for (
int i = 0; i < matrices.size(); i++) {
211 matrices[i].transpose();
217 matrices[i].transpose();
221 void GlobalMasterSymmetry::parseAtoms(
const char *file,
int numTotalAtoms,
int symfileindex) {
222 DebugM(3,
"parseAtoms called\n");
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!");
235 tmdpdb.get_all_positions(atompos);
236 if (K){symmetrykfile = file;}
237 PDB kpdb(symmetrykfile);
241 #ifdef MEM_OPT_VERSION
242 PDBCoreData *atom = tmdpdb.atom(i);
244 PDBAtom *atom = tmdpdb.atom(i);
249 modifyRequestedAtoms().add(i);
251 #ifdef MEM_OPT_VERSION
252 PDBCoreData *atomk = kpdb.atom(i);
260 arr[0] = atompos[i].
x;
261 arr[1] = atompos[i].
y;
262 arr[2] = atompos[i].
z;
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);
275 if (simit != simmap.end()){
276 if (std::find(simit->second.begin(), simit->second.end(), atom->
temperaturefactor()) == simit->second.end()){
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!");
300 void GlobalMasterSymmetry::determineAverage() {
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];
307 delit->second.erase(delit->second.begin(), delit->second.end());
310 map <int, BigReal *>::iterator posit;
311 map <int, vector<int> >::iterator simit = simmap.begin();
312 for (; simit != simmap.end(); ++simit){
315 map <int, BigReal *>::iterator pit = posmap.begin();
316 for (; pit != posmap.end(); ++pit){
delete [] pit->second;}
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()){
325 arr[0] = positions[dit->second[i]].x;
326 arr[1] = positions[dit->second[i]].y;
327 arr[2] = positions[dit->second[i]].z;
329 posmap[dit->second[i]] = arr;
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();
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];
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);
364 void GlobalMasterSymmetry::backTransform(){
365 map <int, BigReal *>::iterator bit = backavg.begin();
366 for (; bit != backavg.end(); ++bit){
delete [] bit->second;}
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();
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];
382 matrices[it->first-1].transpose();
383 gluInvertMatrix(matrices[it->first-1].mat, inverse);
386 matrices[it->first-1].transpose();
389 inv.multpoint(avg+3*k);
391 backavg[it->first] = avg;
396 map <int, BigReal *>::iterator pit = posmap.begin();
397 for (; pit != posmap.end(); ++pit){
delete pit->second;}
399 void GlobalMasterSymmetry::calculate() {
401 modifyAppliedForces().resize(0);
402 modifyForcedAtoms().resize(0);
405 if (currentStep < firstStep || (currentStep >= lastStep && lastStep != -1)) {
417 for ( ; a_i != a_e; ++a_i, ++p_i ){
418 positions[*a_i] = *p_i;
421 map <int, vector<int> >::iterator it;
445 for (it = dmap.begin(); it != dmap.end(); ++it){
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;
456 BigReal *tmpavg = backavg[it->first];
458 for (
int i=0; i<it->second.size(); i++) {
462 k = kdmap[it->first][it->second[i]];
465 k = K/it->second.size();
471 if (currentStep < firstFullStep){
473 (firstFullStep-firstStep);
474 frac = frac*(linear_evolve);
476 if (currentStep > lastFullStep)
479 (lastStep-lastFullStep);
480 frac = frac*(1-linear_evolve);
484 BigReal dx = curpos[3*i ] - tmpavg[3*i];
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;
const Elem * const_iterator
void NAMD_err(const char *err_msg)
void multmatrix(const Matrix4Symmetry &)
premultiply the matrix by the given matrix, this->other * this
BigReal temperaturefactor(void)
void translate(BigReal x, BigReal y, BigReal z)
SimParameters * simParameters
std::ostream & endi(std::ostream &s)
std::ostream & iWARN(std::ostream &s)
void NAMD_die(const char *err_msg)
BigReal MatrixFitRMS(int n, BigReal *v1, BigReal *v2, const BigReal *wt, BigReal *ttt)
StringList * find(const char *name) const
int symmetryFirstFullStep