OOPSE/libmdtools/Bond.cpp

#include "SRI.hpp"
#include "Atom.hpp"
#include <math.h>
#include <iostream>
#include <stdlib.h>


Bond::Bond(){
  
  c_constraint = NULL;
  c_is_constrained = 0;
}

void Bond::set_atoms( Atom &a, Atom &b ){
  
  c_p_a = &a;
  c_p_b = &b;
}  

void Bond::constrain(double bond_distance){

  double dsqr = bond_distance * bond_distance;
  
  c_is_constrained = 1;
  
  c_constraint = new Constraint();
  c_constraint->set_a( c_p_a->getIndex() );
  c_constraint->set_b( c_p_b->getIndex() );
  c_constraint->set_dsqr( dsqr );
}

Bond::~Bond(){
  delete c_constraint;
  c_constraint = 0;
}

void Bond::calc_forces(){
  
  /* return 0 if the bond is constrained and stop wasting cpu */
  
  if(c_is_constrained){
    
    c_potential_E = 0.0;
    return;
  }
  
  vect r_ab; /*the vector whose origin is a and end is b */
  double force; /* the force scaling factor. */
  double Fab_x; /*the x,y, and z components of the force */
  double Fab_y;
  double Fab_z;

  double aR[3], bR[3];
  double aF[3], bF[3];

  /* initialize the vector */
  
  c_p_a->getPos(aR);
  c_p_b->getPos(bR);

  r_ab.x = bR[0] - aR[0];
  r_ab.y = bR[1] - aR[1];
  r_ab.z = bR[2] - aR[2];

  r_ab.length = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
  
  /* calculate the force here */

  force = bond_force(r_ab.length);
  
  Fab_x = -force *  r_ab.x / r_ab.length;
  Fab_y = -force *  r_ab.y / r_ab.length;
  Fab_z = -force *  r_ab.z / r_ab.length;

  aF[0] = Fab_x;
  aF[1] = Fab_y;
  aF[2] = Fab_z;

  bF[0] = -Fab_x;
  bF[1] = -Fab_y;
  bF[2] = -Fab_z;

  c_p_a->addFrc(aF);
  c_p_b->addFrc(bF);

  return;
}
Revision:	670
Committed:	Thu Aug 7 21:47:18 2003 UTC (22 years, 2 months ago) by mmeineke
File size:	1602 byte(s)
Log Message:	switched SimInfo to use a system configuration from SimState rather than arrays from Atom
#	User	Rev	Content
1	mmeineke	377	#include "SRI.hpp"
2			#include "Atom.hpp"
3			#include <math.h>
4			#include <iostream>
5			#include <stdlib.h>
6
7
8
9			Bond::Bond(){
10
11			c_constraint = NULL;
12			c_is_constrained = 0;
13			}
14
15			void Bond::set_atoms( Atom &a, Atom &b ){
16
17			c_p_a = &a;
18			c_p_b = &b;
19			}
20
21			void Bond::constrain(double bond_distance){
22
23			double dsqr = bond_distance * bond_distance;
24
25			c_is_constrained = 1;
26
27			c_constraint = new Constraint();
28			c_constraint->set_a( c_p_a->getIndex() );
29			c_constraint->set_b( c_p_b->getIndex() );
30			c_constraint->set_dsqr( dsqr );
31			}
32
33			Bond::~Bond(){
34			delete c_constraint;
35			c_constraint = 0;
36			}
37
38			void Bond::calc_forces(){
39
40			/* return 0 if the bond is constrained and stop wasting cpu */
41
42			if(c_is_constrained){
43
44			c_potential_E = 0.0;
45			return;
46			}
47
48			vect r_ab; /the vector whose origin is a and end is b /
49			double force; /* the force scaling factor. */
50			double Fab_x; /the x,y, and z components of the force /
51			double Fab_y;
52			double Fab_z;
53
54	mmeineke	670	double aR[3], bR[3];
55			double aF[3], bF[3];
56
57	mmeineke	377	/* initialize the vector */
58	mmeineke	670
59			c_p_a->getPos(aR);
60			c_p_b->getPos(bR);
61	mmeineke	377
62	mmeineke	670	r_ab.x = bR[0] - aR[0];
63			r_ab.y = bR[1] - aR[1];
64			r_ab.z = bR[2] - aR[2];
65
66	mmeineke	377	r_ab.length = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
67
68			/* calculate the force here */
69
70			force = bond_force(r_ab.length);
71
72			Fab_x = -force * r_ab.x / r_ab.length;
73			Fab_y = -force * r_ab.y / r_ab.length;
74			Fab_z = -force * r_ab.z / r_ab.length;
75
76	mmeineke	670	aF[0] = Fab_x;
77			aF[1] = Fab_y;
78			aF[2] = Fab_z;
79	mmeineke	377
80	mmeineke	670	bF[0] = -Fab_x;
81			bF[1] = -Fab_y;
82			bF[2] = -Fab_z;
83	mmeineke	377
84	mmeineke	670	c_p_a->addFrc(aF);
85			c_p_b->addFrc(bF);
86
87	mmeineke	377	return;
88			}