src/restraints/MolecularRestraint.cpp

 /*
 * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#include "restraints/MolecularRestraint.hpp"
#include "math/SquareMatrix3.hpp"
#include "math/SVD.hpp"
#include <utility>

//using namespace JAMA;

namespace oopse {

  void MolecularRestraint::calcForce(std::vector<Vector3d> struc, 
                                     Vector3d molCom){

    assert(struc.size() == ref_.size());

    std::vector<Vector3d>::iterator it;

    // clear out initial values:
    pot_ = 0.0;
    for(it = forces_.begin(); it != forces_.end(); ++it)
      (*it) = 0.0;

   
    if (restType_ & rtDisplacement) {
      Vector3d del = molCom - refCom_;     
      
      RealType r = del.length();
      RealType p = 0.5 * kDisp_ * r * r;

      pot_ += p;

      restInfo_[rtDisplacement] = std::make_pair(r, p);

      for(it = forces_.begin(); it != forces_.end(); ++it)
        (*it) = -kDisp_ * del * scaleFactor_;
    }

    for(it = struc.begin(); it != struc.end(); ++it)
      (*it) -= molCom;

    // rtDisplacement = 1, so anything higher than that requires orientations:
    if (restType_ > 1) {
      Vector3d tBody(0.0);
      
      Mat3x3d R(0.0);
      
      for (int n = 0; n < struc.size(); n++){

        /*
         * correlation matrix R:   
         *   R(i,j) = sum(over n): y(n,i) * x(n,j)  
         *   where x(n) and y(n) are two vector sets   
         */
        
        R += outProduct(struc[n], ref_[n]);
      }    

      // SVD class uses dynamic matrices, so we must wrap the correlation
      // matrix before calling SVD and then unwrap the results into Mat3x3d 
      // and Vector3d before we use them.
      
      DynamicRectMatrix<RealType> Rtmp(3, 3, 0.0);
      DynamicRectMatrix<RealType> vtmp(3, 3);
      DynamicVector<RealType>     stmp(3);
      DynamicRectMatrix<RealType> wtmp(3, 3);

      Rtmp.setSubMatrix(0, 0, R);
    
      // Heavy lifting goes here:
      
      JAMA::SVD<RealType> svd(Rtmp);

      svd.getU(vtmp);
      svd.getSingularValues(stmp);
      svd.getV(wtmp);

      Mat3x3d v;
      Vector3d s;
      Mat3x3d w_tr;
      
      vtmp.getSubMatrix(0, 0, v);
      stmp.getSubVector(0, s);
      wtmp.getSubMatrix(0, 0, w_tr);

      bool is_reflection = (v.determinant() * w_tr.determinant()) < 0.0;
      
      if (is_reflection){        
        v(2, 0) = -v(2, 0);
        v(2, 1) = -v(2, 1);
        v(2, 2) = -v(2, 2);
      }
      
      RotMat3x3d Atrans = v * w_tr.transpose();
      RotMat3x3d A = Atrans.transpose();

      Vector3d eularAngles = A.toEulerAngles();


      RealType twistAngle, swingAngle;
      Vector3d swingAxis;

      Quat4d quat = A.toQuaternion();  

      quat.getTwistSwingAxisAngle(twistAngle, swingAngle, swingAxis);

      RealType tw, sx, sy, ttw, swingX, swingY;
      quat.toTwistSwing(tw, sx, sy);
      quat.toSwingTwist(swingX, swingY, ttw);

      // std::cerr << eularAngles << "\t[" << twistAngle << "," << swingAngle << 
      // "]\t[" << tw << "," << sx << "," << sy << "]\t[" << ttw <<
      // "," << ssx << "," << ssy << "]" << std::endl;

      RealType dVdtwist, dVdswing, dVdswingX, dVdswingY;
      RealType dTwist, dSwing, dSwingX, dSwingY;
      RealType p;

      if (restType_ & rtTwist){
        dTwist = twistAngle - twist0_;
        dVdtwist = kTwist_ * sin(dTwist) ;
        p = kTwist_ * (1.0 - cos(dTwist) ) ;
        pot_ += p;
        tBody -= dVdtwist * V3Z;
        restInfo_[rtTwist] = std::make_pair(twistAngle, p);
      }

//       if (restType_ & rtSwing){
//         dSwing = swingAngle - swing0_;
//         dVdswing = kSwing_ * 2.0 * sin(2.0 * dSwing);
//         p = kSwing_ * (1.0 - cos(2.0 * dSwing));
//         pot_ += p;
//         tBody -= dVdswing * swingAxis;
//         restInfo_[rtSwing] = std::make_pair(swingAngle, p);
//       }

      if (restType_ & rtSwingX){
        dSwingX = swingX - swingX0_;
        dVdswingX = kSwingX_ * 2.0 * sin(2.0 * dSwingX);
        p = kSwingX_ * (1.0 - cos(2.0 * dSwingX));
        pot_ += p;
        tBody -= dVdswingX * V3X;
        restInfo_[rtSwingX] = std::make_pair(swingX, p);
      }
      if (restType_ & rtSwingY){
        dSwingY = swingY - swingY0_;
        dVdswingY = kSwingY_ * 2.0 * sin(2.0 * dSwingY);
        p = kSwingY_ * (1.0 - cos(2.0 * dSwingY));
        pot_ += p;
        tBody -= dVdswingY * V3Y;
        restInfo_[rtSwingY] = std::make_pair(swingY, p);
      }

            
      RealType t2 = dot(tBody, tBody);     
      
      Vector3d rLab, rBody, txr, fBody, fLab;

      for (int i = 0; i < struc.size(); i++) {
                   
        rLab = struc[i];        
        rBody = A * rLab;
        
        txr = cross(tBody, rBody);
        fBody = txr * t2;
        fLab = Atrans * fBody;        
        fLab *= scaleFactor_;
          
        forces_[i] += fLab;
      }

      // test the force vectors and see if it is the right orientation
//       std::cout << struc.size() << std::endl << std::endl;
//       for (int i = 0; i != struc.size(); ++i){
//         std::cout << "H\t" << struc[i].x() << "\t" << struc[i].y() << "\t" << struc[i].z() << "\t";
//         std::cout << forces_[i].x() << "\t" << forces_[i].y() << "\t" << forces_[i].z() << std::endl;
//       }
    }
  }
}
Revision:	1360
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path:	*trunk/src/restraints/MolecularRestraint.cpp*
File size:	7236 byte(s)
Log Message:	Added new restraint infrastructure Added MolecularRestraints Added ObjectRestraints Added RestraintStamp Updated thermodynamic integration to use ObjectRestraints Added Quaternion mathematics for twist swing decompositions Significantly updated RestWriter and RestReader to use dump-like files Added selections for x, y, and z coordinates of atoms Removed monolithic Restraints class Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody Added some rotational capabilities to prinicpalAxisCalculator
#	User	Rev	Content
1	cli2	1360	/*
2			* Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42			#include "restraints/MolecularRestraint.hpp"
43			#include "math/SquareMatrix3.hpp"
44			#include "math/SVD.hpp"
45			#include <utility>
46
47			//using namespace JAMA;
48
49			namespace oopse {
50
51			void MolecularRestraint::calcForce(std::vector<Vector3d> struc,
52			Vector3d molCom){
53
54			assert(struc.size() == ref_.size());
55
56			std::vector<Vector3d>::iterator it;
57
58			// clear out initial values:
59			pot_ = 0.0;
60			for(it = forces_.begin(); it != forces_.end(); ++it)
61			(*it) = 0.0;
62
63
64			if (restType_ & rtDisplacement) {
65			Vector3d del = molCom - refCom_;
66
67			RealType r = del.length();
68			RealType p = 0.5 * kDisp_ * r * r;
69
70			pot_ += p;
71
72			restInfo_[rtDisplacement] = std::make_pair(r, p);
73
74			for(it = forces_.begin(); it != forces_.end(); ++it)
75			(it) = -kDisp_ del * scaleFactor_;
76			}
77
78			for(it = struc.begin(); it != struc.end(); ++it)
79			(*it) -= molCom;
80
81			// rtDisplacement = 1, so anything higher than that requires orientations:
82			if (restType_ > 1) {
83			Vector3d tBody(0.0);
84
85			Mat3x3d R(0.0);
86
87			for (int n = 0; n < struc.size(); n++){
88
89			/*
90			* correlation matrix R:
91			* R(i,j) = sum(over n): y(n,i) * x(n,j)
92			* where x(n) and y(n) are two vector sets
93			*/
94
95			R += outProduct(struc[n], ref_[n]);
96			}
97
98			// SVD class uses dynamic matrices, so we must wrap the correlation
99			// matrix before calling SVD and then unwrap the results into Mat3x3d
100			// and Vector3d before we use them.
101
102			DynamicRectMatrix<RealType> Rtmp(3, 3, 0.0);
103			DynamicRectMatrix<RealType> vtmp(3, 3);
104			DynamicVector<RealType> stmp(3);
105			DynamicRectMatrix<RealType> wtmp(3, 3);
106
107			Rtmp.setSubMatrix(0, 0, R);
108
109			// Heavy lifting goes here:
110
111			JAMA::SVD<RealType> svd(Rtmp);
112
113			svd.getU(vtmp);
114			svd.getSingularValues(stmp);
115			svd.getV(wtmp);
116
117			Mat3x3d v;
118			Vector3d s;
119			Mat3x3d w_tr;
120
121			vtmp.getSubMatrix(0, 0, v);
122			stmp.getSubVector(0, s);
123			wtmp.getSubMatrix(0, 0, w_tr);
124
125			bool is_reflection = (v.determinant() * w_tr.determinant()) < 0.0;
126
127			if (is_reflection){
128			v(2, 0) = -v(2, 0);
129			v(2, 1) = -v(2, 1);
130			v(2, 2) = -v(2, 2);
131			}
132
133			RotMat3x3d Atrans = v * w_tr.transpose();
134			RotMat3x3d A = Atrans.transpose();
135
136			Vector3d eularAngles = A.toEulerAngles();
137
138
139			RealType twistAngle, swingAngle;
140			Vector3d swingAxis;
141
142			Quat4d quat = A.toQuaternion();
143
144			quat.getTwistSwingAxisAngle(twistAngle, swingAngle, swingAxis);
145
146			RealType tw, sx, sy, ttw, swingX, swingY;
147			quat.toTwistSwing(tw, sx, sy);
148			quat.toSwingTwist(swingX, swingY, ttw);
149
150			// std::cerr << eularAngles << "\t[" << twistAngle << "," << swingAngle <<
151			// "]\t[" << tw << "," << sx << "," << sy << "]\t[" << ttw <<
152			// "," << ssx << "," << ssy << "]" << std::endl;
153
154			RealType dVdtwist, dVdswing, dVdswingX, dVdswingY;
155			RealType dTwist, dSwing, dSwingX, dSwingY;
156			RealType p;
157
158			if (restType_ & rtTwist){
159			dTwist = twistAngle - twist0_;
160			dVdtwist = kTwist_ * sin(dTwist) ;
161			p = kTwist_ * (1.0 - cos(dTwist) ) ;
162			pot_ += p;
163			tBody -= dVdtwist * V3Z;
164			restInfo_[rtTwist] = std::make_pair(twistAngle, p);
165			}
166
167			// if (restType_ & rtSwing){
168			// dSwing = swingAngle - swing0_;
169			// dVdswing = kSwing_ * 2.0 * sin(2.0 * dSwing);
170			// p = kSwing_ * (1.0 - cos(2.0 * dSwing));
171			// pot_ += p;
172			// tBody -= dVdswing * swingAxis;
173			// restInfo_[rtSwing] = std::make_pair(swingAngle, p);
174			// }
175
176			if (restType_ & rtSwingX){
177			dSwingX = swingX - swingX0_;
178			dVdswingX = kSwingX_ * 2.0 * sin(2.0 * dSwingX);
179			p = kSwingX_ * (1.0 - cos(2.0 * dSwingX));
180			pot_ += p;
181			tBody -= dVdswingX * V3X;
182			restInfo_[rtSwingX] = std::make_pair(swingX, p);
183			}
184			if (restType_ & rtSwingY){
185			dSwingY = swingY - swingY0_;
186			dVdswingY = kSwingY_ * 2.0 * sin(2.0 * dSwingY);
187			p = kSwingY_ * (1.0 - cos(2.0 * dSwingY));
188			pot_ += p;
189			tBody -= dVdswingY * V3Y;
190			restInfo_[rtSwingY] = std::make_pair(swingY, p);
191			}
192
193
194			RealType t2 = dot(tBody, tBody);
195
196			Vector3d rLab, rBody, txr, fBody, fLab;
197
198			for (int i = 0; i < struc.size(); i++) {
199
200			rLab = struc[i];
201			rBody = A * rLab;
202
203			txr = cross(tBody, rBody);
204			fBody = txr * t2;
205			fLab = Atrans * fBody;
206			fLab *= scaleFactor_;
207
208			forces_[i] += fLab;
209			}
210
211			// test the force vectors and see if it is the right orientation
212			// std::cout << struc.size() << std::endl << std::endl;
213			// for (int i = 0; i != struc.size(); ++i){
214			// std::cout << "H\t" << struc[i].x() << "\t" << struc[i].y() << "\t" << struc[i].z() << "\t";
215			// std::cout << forces_[i].x() << "\t" << forces_[i].y() << "\t" << forces_[i].z() << std::endl;
216			// }
217			}
218			}
219			}