src/primitives/DirectionalAtom.cpp

/*
 * Copyright (c) 2005 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 "primitives/DirectionalAtom.hpp"
#include "utils/simError.h"
namespace oopse {
  
  DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType) 
    : Atom(dAtomType){
    objType_= otDAtom;
    if (dAtomType->isMultipole()) {
      electroBodyFrame_ = dAtomType->getElectroBodyFrame();
    }
    
    // Check if one of the diagonal inertia tensor of this directional
    // atom is zero:
    int nLinearAxis = 0;
    Mat3x3d inertiaTensor = getI();
    for (int i = 0; i < 3; i++) {    
      if (fabs(inertiaTensor(i, i)) < oopse::epsilon) {
        linear_ = true;
        linearAxis_ = i;
        ++ nLinearAxis;
      }
    }

    if (nLinearAxis > 1) {
      sprintf( painCave.errMsg,
               "Directional Atom warning.\n"
               "\tOOPSE found more than one axis in this directional atom with a vanishing \n"
               "\tmoment of inertia.");
      painCave.isFatal = 0;
      simError();
    }    
  }
  
  Mat3x3d DirectionalAtom::getI() {
    return static_cast<DirectionalAtomType*>(getAtomType())->getI();
  }    
  
  void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
    ((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
    if (atomType_->isMultipole()) {
      ((snapshotMan_->getPrevSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
    }
  }
  
  
  void DirectionalAtom::setA(const RotMat3x3d& a) {
    ((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
    
    if (atomType_->isMultipole()) {
      ((snapshotMan_->getCurrentSnapshot())->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;
    }
  }    
  
  void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
    ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
    
    if (atomType_->isMultipole()) {
      ((snapshotMan_->getSnapshot(snapshotNo))->*storage_).electroFrame[localIndex_] = a.transpose() * electroBodyFrame_;    
    }
  }    
  
  void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
    setA(m *getA());
  }
  
  std::vector<RealType> DirectionalAtom::getGrad() {
    std::vector<RealType> grad(6, 0.0);
    Vector3d force;
    Vector3d torque;
    Vector3d myEuler;
    RealType phi, theta, psi;
    RealType cphi, sphi, ctheta, stheta;
    Vector3d ephi;
    Vector3d etheta;
    Vector3d epsi;
    
    force = getFrc();
    torque =getTrq();
    myEuler = getA().toEulerAngles();
    
    phi = myEuler[0];
    theta = myEuler[1];
    psi = myEuler[2];
    
    cphi = cos(phi);
    sphi = sin(phi);
    ctheta = cos(theta);
    stheta = sin(theta);
    
    // get unit vectors along the phi, theta and psi rotation axes
    
    ephi[0] = 0.0;
    ephi[1] = 0.0;
    ephi[2] = 1.0;
    
    etheta[0] = -sphi;
    etheta[1] =  cphi;
    etheta[2] =  0.0;
    
    epsi[0] = stheta * cphi;
    epsi[1] = stheta * sphi;
    epsi[2] = ctheta;
    
    //gradient is equal to -force
    for (int j = 0 ; j<3; j++)
      grad[j] = -force[j];
    
    for (int j = 0; j < 3; j++ ) {      
      grad[3] -= torque[j]*ephi[j];
      grad[4] -= torque[j]*etheta[j];
      grad[5] -= torque[j]*epsi[j];      
    }
    
    return grad;
  }    
  
  void DirectionalAtom::accept(BaseVisitor* v) {
    v->visit(this);
  }
}

Revision:	1360
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path:	*trunk/src/primitives/DirectionalAtom.cpp*
File size:	5370 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	gezelter	507	/*
2	gezelter	246	* Copyright (c) 2005 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	tim	3	#include "primitives/DirectionalAtom.hpp"
43	tim	273	#include "utils/simError.h"
44	gezelter	246	namespace oopse {
45	gezelter	1211
46	gezelter	507	DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47			: Atom(dAtomType){
48	gezelter	1211	objType_= otDAtom;
49			if (dAtomType->isMultipole()) {
50			electroBodyFrame_ = dAtomType->getElectroBodyFrame();
51			}
52
53			// Check if one of the diagonal inertia tensor of this directional
54			// atom is zero:
55			int nLinearAxis = 0;
56			Mat3x3d inertiaTensor = getI();
57			for (int i = 0; i < 3; i++) {
58			if (fabs(inertiaTensor(i, i)) < oopse::epsilon) {
59			linear_ = true;
60			linearAxis_ = i;
61			++ nLinearAxis;
62	gezelter	507	}
63	tim	273	}
64	gezelter	2
65	gezelter	1211	if (nLinearAxis > 1) {
66			sprintf( painCave.errMsg,
67			"Directional Atom warning.\n"
68			"\tOOPSE found more than one axis in this directional atom with a vanishing \n"
69			"\tmoment of inertia.");
70			painCave.isFatal = 0;
71			simError();
72			}
73			}
74
75	gezelter	507	Mat3x3d DirectionalAtom::getI() {
76	gezelter	246	return static_cast<DirectionalAtomType*>(getAtomType())->getI();
77	gezelter	507	}
78	gezelter	1211
79	gezelter	507	void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
80	gezelter	246	((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
81			if (atomType_->isMultipole()) {
82	gezelter	507	((snapshotMan_->getPrevSnapshot())->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
83	gezelter	205	}
84	gezelter	507	}
85	gezelter	1211
86
87	gezelter	507	void DirectionalAtom::setA(const RotMat3x3d& a) {
88	gezelter	246	((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
89	gezelter	1211
90	gezelter	246	if (atomType_->isMultipole()) {
91	gezelter	507	((snapshotMan_->getCurrentSnapshot())->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
92	gezelter	2	}
93	gezelter	507	}
94	gezelter	1211
95	gezelter	507	void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
96	gezelter	246	((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
97	gezelter	1211
98	gezelter	246	if (atomType_->isMultipole()) {
99	gezelter	507	((snapshotMan_->getSnapshot(snapshotNo))->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
100	gezelter	2	}
101	gezelter	507	}
102	gezelter	1211
103	gezelter	507	void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
104	gezelter	246	setA(m *getA());
105	gezelter	507	}
106	gezelter	1211
107	tim	963	std::vector<RealType> DirectionalAtom::getGrad() {
108			std::vector<RealType> grad(6, 0.0);
109	gezelter	246	Vector3d force;
110			Vector3d torque;
111			Vector3d myEuler;
112	tim	963	RealType phi, theta, psi;
113			RealType cphi, sphi, ctheta, stheta;
114	gezelter	246	Vector3d ephi;
115			Vector3d etheta;
116			Vector3d epsi;
117	gezelter	1211
118	gezelter	246	force = getFrc();
119			torque =getTrq();
120			myEuler = getA().toEulerAngles();
121	gezelter	1211
122	gezelter	246	phi = myEuler[0];
123			theta = myEuler[1];
124			psi = myEuler[2];
125	gezelter	1211
126	gezelter	246	cphi = cos(phi);
127			sphi = sin(phi);
128			ctheta = cos(theta);
129			stheta = sin(theta);
130	gezelter	1211
131	gezelter	246	// get unit vectors along the phi, theta and psi rotation axes
132	gezelter	1211
133	gezelter	246	ephi[0] = 0.0;
134			ephi[1] = 0.0;
135			ephi[2] = 1.0;
136	gezelter	1211
137	cli2	1360	etheta[0] = -sphi;
138			etheta[1] = cphi;
139			etheta[2] = 0.0;
140	gezelter	1211
141	gezelter	246	epsi[0] = stheta * cphi;
142			epsi[1] = stheta * sphi;
143			epsi[2] = ctheta;
144	gezelter	1211
145	gezelter	246	//gradient is equal to -force
146			for (int j = 0 ; j<3; j++)
147	gezelter	507	grad[j] = -force[j];
148	gezelter	1211
149			for (int j = 0; j < 3; j++ ) {
150	tim	642	grad[3] -= torque[j]*ephi[j];
151			grad[4] -= torque[j]*etheta[j];
152	gezelter	1211	grad[5] -= torque[j]*epsi[j];
153	gezelter	246	}
154	gezelter	2
155	gezelter	246	return grad;
156	gezelter	507	}
157	gezelter	1211
158	gezelter	507	void DirectionalAtom::accept(BaseVisitor* v) {
159	gezelter	246	v->visit(this);
160	gezelter	1211	}
161	gezelter	2	}
162