OOPSE/libmdtools/Shake.cpp

#include "Shake.hpp"
#include "SimInfo.hpp"
#include <cmath>
#include "simError.h"
////////////////////////////////////////////////////////////////////////////////
//Implementation of DCShakeFunctor
////////////////////////////////////////////////////////////////////////////////
int DCShakeFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
  double posA[3];
  double posB[3];
  double oldPosA[3];
  double oldPosB[3];
  double velA[3];
  double velB[3];
  double pab[3];
  double rab[3];
  double rma, rmb;
  double dx, dy, dz;
  double rpab;
  double rabsq, pabsq, rpabsq;
  double diffsq;
  double gab;
  double dt;
  
  dt = info->dt;
  
  consAtom1->getPos(posA);
  consAtom2->getPos(posB);

   
  pab[0] = posA[0] - posB[0];
  pab[1] = posA[1] - posB[1];
  pab[2] = posA[2] - posB[2];

  //periodic boundary condition

  info->wrapVector(pab);

  pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];

  rabsq = curPair->getBondLength2();
  diffsq = rabsq - pabsq;

  // the original rattle code from alan tidesley
  if (fabs(diffsq) > (consTolerance * rabsq * 2)){
    
    consAtom1->getOldPos(oldPosA);
    consAtom2->getOldPos(oldPosB);      
    
    rab[0] = oldPosA[0] - oldPosB[0];
    rab[1] = oldPosA[1] - oldPosB[1];
    rab[2] = oldPosA[2] - oldPosB[2];

    info->wrapVector(rab);

    rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];

    rpabsq = rpab * rpab;


    if (rpabsq < (rabsq * -diffsq)){
      return consFail;
    }

    rma = 1.0 / consAtom1->getMass();
    rmb = 1.0 / consAtom2->getMass();

    gab = diffsq / (2.0 * (rma + rmb) * rpab);

    dx = rab[0] * gab;
    dy = rab[1] * gab;
    dz = rab[2] * gab;

    posA[0] += rma * dx;
    posA[1] += rma * dy;
    posA[2] += rma * dz;

    consAtom1->setPos(posA);

    posB[0] -= rmb * dx;
    posB[1] -= rmb * dy;
    posB[2] -= rmb * dz;

    consAtom2->setPos(posB);

    dx = dx / dt;
    dy = dy / dt;
    dz = dz / dt;

    consAtom1->getVel(velA);

    velA[0] += rma * dx;
    velA[1] += rma * dy;
    velA[2] += rma * dz;

    consAtom1->setVel(velA);

    consAtom2->getVel(velB);

    velB[0] -= rmb * dx;
    velB[1] -= rmb * dy;
    velB[2] -= rmb * dz;

    consAtom2->setVel(velB);

    return consSuccess;
  }
  else
    return consAlready;
  
}


int DCShakeFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
  return consElemHandlerFail;
} 

int DCShakeFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
  return consElemHandlerFail;
}

////////////////////////////////////////////////////////////////////////////////
//Implementation of JCShakeFunctor
////////////////////////////////////////////////////////////////////////////////
int JCShakeFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
  return consElemHandlerFail;
}


int JCShakeFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
  return consElemHandlerFail;

}
 
int JCShakeFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
  return consElemHandlerFail;
}
Revision:	1232
Committed:	Thu Jun 3 21:51:55 2004 UTC (21 years, 5 months ago) by tim
File size:	3100 byte(s)
Log Message:	new implementation of constraint
#	User	Rev	Content
1	tim	1232	#include "Shake.hpp"
2			#include "SimInfo.hpp"
3			#include <cmath>
4			#include "simError.h"
5			////////////////////////////////////////////////////////////////////////////////
6			//Implementation of DCShakeFunctor
7			////////////////////////////////////////////////////////////////////////////////
8			int DCShakeFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
9			double posA[3];
10			double posB[3];
11			double oldPosA[3];
12			double oldPosB[3];
13			double velA[3];
14			double velB[3];
15			double pab[3];
16			double rab[3];
17			double rma, rmb;
18			double dx, dy, dz;
19			double rpab;
20			double rabsq, pabsq, rpabsq;
21			double diffsq;
22			double gab;
23			double dt;
24
25			dt = info->dt;
26
27			consAtom1->getPos(posA);
28			consAtom2->getPos(posB);
29
30
31			pab[0] = posA[0] - posB[0];
32			pab[1] = posA[1] - posB[1];
33			pab[2] = posA[2] - posB[2];
34
35			//periodic boundary condition
36
37			info->wrapVector(pab);
38
39			pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
40
41			rabsq = curPair->getBondLength2();
42			diffsq = rabsq - pabsq;
43
44			// the original rattle code from alan tidesley
45			if (fabs(diffsq) > (consTolerance * rabsq * 2)){
46
47			consAtom1->getOldPos(oldPosA);
48			consAtom2->getOldPos(oldPosB);
49
50			rab[0] = oldPosA[0] - oldPosB[0];
51			rab[1] = oldPosA[1] - oldPosB[1];
52			rab[2] = oldPosA[2] - oldPosB[2];
53
54			info->wrapVector(rab);
55
56			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
57
58			rpabsq = rpab * rpab;
59
60
61			if (rpabsq < (rabsq * -diffsq)){
62			return consFail;
63			}
64
65			rma = 1.0 / consAtom1->getMass();
66			rmb = 1.0 / consAtom2->getMass();
67
68			gab = diffsq / (2.0 * (rma + rmb) * rpab);
69
70			dx = rab[0] * gab;
71			dy = rab[1] * gab;
72			dz = rab[2] * gab;
73
74			posA[0] += rma * dx;
75			posA[1] += rma * dy;
76			posA[2] += rma * dz;
77
78			consAtom1->setPos(posA);
79
80			posB[0] -= rmb * dx;
81			posB[1] -= rmb * dy;
82			posB[2] -= rmb * dz;
83
84			consAtom2->setPos(posB);
85
86			dx = dx / dt;
87			dy = dy / dt;
88			dz = dz / dt;
89
90			consAtom1->getVel(velA);
91
92			velA[0] += rma * dx;
93			velA[1] += rma * dy;
94			velA[2] += rma * dz;
95
96			consAtom1->setVel(velA);
97
98			consAtom2->getVel(velB);
99
100			velB[0] -= rmb * dx;
101			velB[1] -= rmb * dy;
102			velB[2] -= rmb * dz;
103
104			consAtom2->setVel(velB);
105
106			return consSuccess;
107			}
108			else
109			return consAlready;
110
111			}
112
113
114			int DCShakeFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
115			return consElemHandlerFail;
116			}
117
118			int DCShakeFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
119			return consElemHandlerFail;
120			}
121
122			////////////////////////////////////////////////////////////////////////////////
123			//Implementation of JCShakeFunctor
124			////////////////////////////////////////////////////////////////////////////////
125			int JCShakeFunctor::operator()(ConstraintAtom* consAtom1, ConstraintAtom* consAtom2){
126			return consElemHandlerFail;
127			}
128
129
130
131			int JCShakeFunctor::operator()(ConstraintAtom* consAtom,ConstraintRigidBody* consRB){
132			return consElemHandlerFail;
133
134			}
135
136			int JCShakeFunctor::operator()(ConstraintRigidBody* consRB1, ConstraintRigidBody* consRB2){
137			return consElemHandlerFail;
138			}