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:	3520
Committed:	Mon Sep 7 16:31:51 2009 UTC (15 years, 10 months ago) by cli2
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	2204	/*
2	gezelter	1930	* 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	1492	#include "primitives/DirectionalAtom.hpp"
43	tim	1957	#include "utils/simError.h"
44	gezelter	1930	namespace oopse {
45	gezelter	3320
46	gezelter	2204	DirectionalAtom::DirectionalAtom(DirectionalAtomType* dAtomType)
47			: Atom(dAtomType){
48	gezelter	3320	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	2204	}
63	tim	1957	}
64	gezelter	1490
65	gezelter	3320	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	2204	Mat3x3d DirectionalAtom::getI() {
76	gezelter	1930	return static_cast<DirectionalAtomType*>(getAtomType())->getI();
77	gezelter	2204	}
78	gezelter	3320
79	gezelter	2204	void DirectionalAtom::setPrevA(const RotMat3x3d& a) {
80	gezelter	1930	((snapshotMan_->getPrevSnapshot())->*storage_).aMat[localIndex_] = a;
81			if (atomType_->isMultipole()) {
82	gezelter	2204	((snapshotMan_->getPrevSnapshot())->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
83	gezelter	1709	}
84	gezelter	2204	}
85	gezelter	3320
86
87	gezelter	2204	void DirectionalAtom::setA(const RotMat3x3d& a) {
88	gezelter	1930	((snapshotMan_->getCurrentSnapshot())->*storage_).aMat[localIndex_] = a;
89	gezelter	3320
90	gezelter	1930	if (atomType_->isMultipole()) {
91	gezelter	2204	((snapshotMan_->getCurrentSnapshot())->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
92	gezelter	1490	}
93	gezelter	2204	}
94	gezelter	3320
95	gezelter	2204	void DirectionalAtom::setA(const RotMat3x3d& a, int snapshotNo) {
96	gezelter	1930	((snapshotMan_->getSnapshot(snapshotNo))->*storage_).aMat[localIndex_] = a;
97	gezelter	3320
98	gezelter	1930	if (atomType_->isMultipole()) {
99	gezelter	2204	((snapshotMan_->getSnapshot(snapshotNo))->storage_).electroFrame[localIndex_] = a.transpose() electroBodyFrame_;
100	gezelter	1490	}
101	gezelter	2204	}
102	gezelter	3320
103	gezelter	2204	void DirectionalAtom::rotateBy(const RotMat3x3d& m) {
104	gezelter	1930	setA(m *getA());
105	gezelter	2204	}
106	gezelter	3320
107	tim	2759	std::vector<RealType> DirectionalAtom::getGrad() {
108			std::vector<RealType> grad(6, 0.0);
109	gezelter	1930	Vector3d force;
110			Vector3d torque;
111			Vector3d myEuler;
112	tim	2759	RealType phi, theta, psi;
113			RealType cphi, sphi, ctheta, stheta;
114	gezelter	1930	Vector3d ephi;
115			Vector3d etheta;
116			Vector3d epsi;
117	gezelter	3320
118	gezelter	1930	force = getFrc();
119			torque =getTrq();
120			myEuler = getA().toEulerAngles();
121	gezelter	3320
122	gezelter	1930	phi = myEuler[0];
123			theta = myEuler[1];
124			psi = myEuler[2];
125	gezelter	3320
126	gezelter	1930	cphi = cos(phi);
127			sphi = sin(phi);
128			ctheta = cos(theta);
129			stheta = sin(theta);
130	gezelter	3320
131	gezelter	1930	// get unit vectors along the phi, theta and psi rotation axes
132	gezelter	3320
133	gezelter	1930	ephi[0] = 0.0;
134			ephi[1] = 0.0;
135			ephi[2] = 1.0;
136	gezelter	3320
137	cli2	3520	etheta[0] = -sphi;
138			etheta[1] = cphi;
139			etheta[2] = 0.0;
140	gezelter	3320
141	gezelter	1930	epsi[0] = stheta * cphi;
142			epsi[1] = stheta * sphi;
143			epsi[2] = ctheta;
144	gezelter	3320
145	gezelter	1930	//gradient is equal to -force
146			for (int j = 0 ; j<3; j++)
147	gezelter	2204	grad[j] = -force[j];
148	gezelter	3320
149			for (int j = 0; j < 3; j++ ) {
150	tim	2341	grad[3] -= torque[j]*ephi[j];
151			grad[4] -= torque[j]*etheta[j];
152	gezelter	3320	grad[5] -= torque[j]*epsi[j];
153	gezelter	1930	}
154	gezelter	1490
155	gezelter	1930	return grad;
156	gezelter	2204	}
157	gezelter	3320
158	gezelter	2204	void DirectionalAtom::accept(BaseVisitor* v) {
159	gezelter	1930	v->visit(this);
160	gezelter	3320	}
161	gezelter	1490	}
162