OOPSE/libmdtools/GhostBend.cpp


#include <math.h>
#include <iostream>
#include <stdlib.h>

#include "simError.h"
#include "SRI.hpp"
#include "Atom.hpp"


GhostBend::GhostBend( Atom &a, Atom &b ){

  c_p_a = &a;
  
  if( !b.isDirectional() ){
    
    // if atom b is not directional, then bad things will happen
    
    sprintf( painCave.errMsg,
             " Ghost Bend error: Atom # %d of type \"%s\" is not "
             "directional.\n",
             b.getIndex(),
             b.getType() );
    painCave.isFatal = 1;
    simError();
  }    

  atomB = ( DirectionalAtom* ) &b;
  
  c_potential_E = 0.0;

#ifdef IS_MPI
  sprintf( checkPointMsg,
           "Succesfull creation of Ghost bend\n" );
  MPIcheckPoint();
#endif // is_mpi

}


void GhostBend::calc_forces(){
  
  double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2;
  double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl;
  
  double sina2, sinai;

  double comf2, comf3, comf4;
  double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz;
  // double dcsjdx, dcsjdy, dcsjdz;
  double dadxi, dadyi, dadzi;
  double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj;
  double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj;
  double u[3];
  
  dx = c_p_a->getX() - atomB->getX();
  dy = c_p_a->getY() - atomB->getY();
  dz = c_p_a->getZ() - atomB->getZ();
 
  atomB->getU(u);

  gx = u[0];
  gy = u[1];
  gz = u[2];
  
  dx2 = dx * dx;
  dy2 = dy * dy;
  dz2 = dz * dz;

  gx2 = gx * gx;
  gy2 = gy * gy;
  gz2 = gz * gz;
  
  rij2 = dx2 + dy2 + dz2;
  rkj2 = gx2 + gy2 + gz2;
  
  riji2 = 1.0 / rij2;
  rkji2 = 1.0 / rkj2;

  dot = dx * gx + dy * gy + dz * gz;
  denom = sqrt((riji2 * rkji2));
  cosang = dot * denom;

  if(cosang > 1.0)cosang = 1.0;
  if(cosang < -1.0) cosang = -1.0;

  angl = acos(cosang);
  angl = angl * 180.0 / M_PI;

  sina2 = 1.0 - cosang*cosang;
  if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12;
  sinai = 1.0 / sqrt(sina2);

  comf2 = cosang * riji2;
  comf3 = cosang * rkji2;
  comf4 = bend_force(angl);

  dcsidx = gx*denom - comf2*dx;
  dcsidy = gy*denom - comf2*dy;
  dcsidz = gz*denom - comf2*dz;
  
  dcskdx = dx*denom - comf3*gx;
  dcskdy = dy*denom - comf3*gy;
  dcskdz = dz*denom - comf3*gz;
  
//   dcsjdx = -dcsidx - dcskdx;
//   dcsjdy = -dcsidy - dcskdy;
//   dcsjdz = -dcsidz - dcskdz;

  dadxi = -sinai*dcsidx;
  dadyi = -sinai*dcsidy;
  dadzi = -sinai*dcsidz;

  dadxk = -sinai*dcskdx;
  dadyk = -sinai*dcskdy;
  dadzk = -sinai*dcskdz;

//   dadxj = -dadxi - dadxk;
//   dadyj = -dadyi - dadyk;
//   dadzj = -dadzi - dadzk;

  daxi = comf4*dadxi;
  dayi = comf4*dadyi;
  dazi = comf4*dadzi;

  daxk = comf4*dadxk;
  dayk = comf4*dadyk;
  dazk = comf4*dadzk;
  
  daxj = -daxi - daxk;
  dayj = -dayi - dayk;
  dazj = -dazi - dazk;

  c_p_a->addFx(daxi);
  c_p_a->addFy(dayi);
  c_p_a->addFz(dazi);

  atomB->addFx(daxj + daxk);
  atomB->addFy(dayj + dayk);
  atomB->addFz(dazj + dazk);

  atomB->addTx(gy*dazk - gz*dayk);
  atomB->addTy(gz*daxk - gx*dazk);
  atomB->addTz(gx*dayk - gy*daxk);

  return;
}

void GhostBend::setConstants( double the_c1, double the_c2, double the_c3, 
                                  double the_Th0 ){
  c1 = the_c1;
  c2 = the_c2;
  c3 = the_c3;
  theta0 = the_Th0;
}


double GhostBend::bend_force( double theta ){

  double dt, dt2;
  double force;

  dt = ( theta - theta0 ) * M_PI / 180.0;
  dt2 = dt * dt;

  c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3;
  force = -( ( 2.0 * c1 * dt ) + c2 );
  return force;
}
Revision:	378
Committed:	Fri Mar 21 17:42:12 2003 UTC (21 years, 3 months ago) by mmeineke
File size:	3365 byte(s)
Log Message:	This commit was generated by cvs2svn to compensate for changes in r377, which included commits to RCS files with non-trunk default branches.
#	User	Rev	Content
1	mmeineke	377
2			#include <math.h>
3			#include <iostream>
4			#include <stdlib.h>
5
6			#include "simError.h"
7			#include "SRI.hpp"
8			#include "Atom.hpp"
9
10
11
12			GhostBend::GhostBend( Atom &a, Atom &b ){
13
14			c_p_a = &a;
15
16			if( !b.isDirectional() ){
17
18			// if atom b is not directional, then bad things will happen
19
20			sprintf( painCave.errMsg,
21			" Ghost Bend error: Atom # %d of type \"%s\" is not "
22			"directional.\n",
23			b.getIndex(),
24			b.getType() );
25			painCave.isFatal = 1;
26			simError();
27			}
28
29			atomB = ( DirectionalAtom* ) &b;
30
31			c_potential_E = 0.0;
32
33			#ifdef IS_MPI
34			sprintf( checkPointMsg,
35			"Succesfull creation of Ghost bend\n" );
36			MPIcheckPoint();
37			#endif // is_mpi
38
39			}
40
41
42			void GhostBend::calc_forces(){
43
44			double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2;
45			double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl;
46
47			double sina2, sinai;
48
49			double comf2, comf3, comf4;
50			double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz;
51			// double dcsjdx, dcsjdy, dcsjdz;
52			double dadxi, dadyi, dadzi;
53			double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj;
54			double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj;
55			double u[3];
56
57			dx = c_p_a->getX() - atomB->getX();
58			dy = c_p_a->getY() - atomB->getY();
59			dz = c_p_a->getZ() - atomB->getZ();
60
61			atomB->getU(u);
62
63			gx = u[0];
64			gy = u[1];
65			gz = u[2];
66
67			dx2 = dx * dx;
68			dy2 = dy * dy;
69			dz2 = dz * dz;
70
71			gx2 = gx * gx;
72			gy2 = gy * gy;
73			gz2 = gz * gz;
74
75			rij2 = dx2 + dy2 + dz2;
76			rkj2 = gx2 + gy2 + gz2;
77
78			riji2 = 1.0 / rij2;
79			rkji2 = 1.0 / rkj2;
80
81			dot = dx * gx + dy * gy + dz * gz;
82			denom = sqrt((riji2 * rkji2));
83			cosang = dot * denom;
84
85			if(cosang > 1.0)cosang = 1.0;
86			if(cosang < -1.0) cosang = -1.0;
87
88			angl = acos(cosang);
89			angl = angl * 180.0 / M_PI;
90
91			sina2 = 1.0 - cosang*cosang;
92			if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12;
93			sinai = 1.0 / sqrt(sina2);
94
95			comf2 = cosang * riji2;
96			comf3 = cosang * rkji2;
97			comf4 = bend_force(angl);
98
99			dcsidx = gxdenom - comf2dx;
100			dcsidy = gydenom - comf2dy;
101			dcsidz = gzdenom - comf2dz;
102
103			dcskdx = dxdenom - comf3gx;
104			dcskdy = dydenom - comf3gy;
105			dcskdz = dzdenom - comf3gz;
106
107			// dcsjdx = -dcsidx - dcskdx;
108			// dcsjdy = -dcsidy - dcskdy;
109			// dcsjdz = -dcsidz - dcskdz;
110
111			dadxi = -sinai*dcsidx;
112			dadyi = -sinai*dcsidy;
113			dadzi = -sinai*dcsidz;
114
115			dadxk = -sinai*dcskdx;
116			dadyk = -sinai*dcskdy;
117			dadzk = -sinai*dcskdz;
118
119			// dadxj = -dadxi - dadxk;
120			// dadyj = -dadyi - dadyk;
121			// dadzj = -dadzi - dadzk;
122
123			daxi = comf4*dadxi;
124			dayi = comf4*dadyi;
125			dazi = comf4*dadzi;
126
127			daxk = comf4*dadxk;
128			dayk = comf4*dadyk;
129			dazk = comf4*dadzk;
130
131			daxj = -daxi - daxk;
132			dayj = -dayi - dayk;
133			dazj = -dazi - dazk;
134
135			c_p_a->addFx(daxi);
136			c_p_a->addFy(dayi);
137			c_p_a->addFz(dazi);
138
139			atomB->addFx(daxj + daxk);
140			atomB->addFy(dayj + dayk);
141			atomB->addFz(dazj + dazk);
142
143			atomB->addTx(gydazk - gzdayk);
144			atomB->addTy(gzdaxk - gxdazk);
145			atomB->addTz(gxdayk - gydaxk);
146
147			return;
148			}
149
150			void GhostBend::setConstants( double the_c1, double the_c2, double the_c3,
151			double the_Th0 ){
152			c1 = the_c1;
153			c2 = the_c2;
154			c3 = the_c3;
155			theta0 = the_Th0;
156			}
157
158
159			double GhostBend::bend_force( double theta ){
160
161			double dt, dt2;
162			double force;
163
164			dt = ( theta - theta0 ) * M_PI / 180.0;
165			dt2 = dt * dt;
166
167			c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3;
168			force = -( ( 2.0 * c1 * dt ) + c2 );
169			return force;
170			}