File: stretchComponent.C

package info (click to toggle)
ball 1.5.0%2Bgit20180813.37fc53c-3
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 239,848 kB
  • sloc: cpp: 326,149; ansic: 4,208; python: 2,303; yacc: 1,778; lex: 1,099; xml: 958; sh: 322; makefile: 93
file content (80 lines) | stat: -rw-r--r-- 1,870 bytes parent folder | download | duplicates (8)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <BALL/MOLMEC/COMMON/stretchComponent.h>

#include <BALL/MOLMEC/COMMON/forceField.h>

namespace BALL
{

	StretchComponent::StretchComponent()
	{
	}

	StretchComponent::StretchComponent(ForceField& ff)
		: ForceFieldComponent(ff)
	{
	}

	StretchComponent::~StretchComponent()
	{
	}

	double StretchComponent::updateEnergy()
	{
		// initial energy is zero
		energy_ = 0;

		bool use_selection = getForceField()->getUseSelection();

		// iterate over all bonds, sum up the energies
		for (Size i = 0; i < stretch_.size(); i++)
		{
			const Atom* atom1 = stretch_[i].atom1;
			const Atom* atom2 = stretch_[i].atom2;
			if (!use_selection || atom1->isSelected() || atom2->isSelected())
			{
				double distance = atom1->getDistance(*atom2);
				energy_ += stretch_[i].values.k * (distance - stretch_[i].values.r0) * (distance - stretch_[i].values.r0);
			}
		}

		return energy_;
	}

	// calculates and adds its forces to the current forces of the force field
	void StretchComponent::updateForces()
	{
		if (getForceField() == 0)
		{
			return;
		}

		bool use_selection = getForceField()->getUseSelection();

		// iterate over all bonds, update the forces
		for (Size i = 0 ; i < stretch_.size(); i++)
		{
			Atom* atom1 = stretch_[i].atom1;
			Atom* atom2 = stretch_[i].atom2;
			Vector3 direction(atom1->getPosition() - atom2->getPosition());
			double distance = direction.getLength();

			if (distance != 0.0)
			{
				// unit conversion: from kJ/(mol A) -> N
				//   kJ -> J: 1e3
				//   A  -> m: 1e10
				//   J/mol -> J: Avogadro
				direction *= 1e13 / Constants::AVOGADRO * 2 * stretch_[i].values.k * (distance - stretch_[i].values.r0) / distance;

				if (!use_selection || atom1->isSelected())
				{
					atom1->getForce() -= direction;
				}
				if (!use_selection || atom2->isSelected())
				{
					atom2->getForce() += direction;
				}
			}
		}
	}
}