Commit ce7baae5 by Francois Gygi

removed trailing blanks in *.[Ch] *.mk Makefile


git-svn-id: http://qboxcode.org/svn/qb/trunk@515 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 14e9be36
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// AndersonMixer.C // AndersonMixer.C
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: AndersonMixer.C,v 1.4 2007-03-17 01:14:00 fgygi Exp $ // $Id: AndersonMixer.C,v 1.5 2007-10-19 16:24:03 fgygi Exp $
#include "AndersonMixer.h" #include "AndersonMixer.h"
#include "blas.h" #include "blas.h"
...@@ -28,13 +28,13 @@ void AndersonMixer::update(const double* f, double* theta, double* fbar) ...@@ -28,13 +28,13 @@ void AndersonMixer::update(const double* f, double* theta, double* fbar)
// tmp0 = delta_F = f - flast // tmp0 = delta_F = f - flast
for ( int i = 0; i < n_; i++ ) for ( int i = 0; i < n_; i++ )
tmp0[i] = f[i] - flast_[i]; tmp0[i] = f[i] - flast_[i];
// a = delta_F * F // a = delta_F * F
double a = ddot(&n_, &tmp0[0], &ione, f, &ione); double a = ddot(&n_, &tmp0[0], &ione, f, &ione);
// b = delta_F * delta_F // b = delta_F * delta_F
double b = ddot(&n_, &tmp0[0], &ione, &tmp0[0], &ione); double b = ddot(&n_, &tmp0[0], &ione, &tmp0[0], &ione);
if ( pctxt_ != 0 ) if ( pctxt_ != 0 )
{ {
double work[2] = { a, b }; double work[2] = { a, b };
...@@ -47,16 +47,16 @@ void AndersonMixer::update(const double* f, double* theta, double* fbar) ...@@ -47,16 +47,16 @@ void AndersonMixer::update(const double* f, double* theta, double* fbar)
*theta = - a / b; *theta = - a / b;
else else
*theta = 0.0; *theta = 0.0;
// test if residual has increased // test if residual has increased
if ( *theta <= -1.0 ) if ( *theta <= -1.0 )
{ {
*theta = theta_nc_; *theta = theta_nc_;
} }
*theta = min(theta_max_,*theta); *theta = min(theta_max_,*theta);
} }
// fbar = f + theta * ( f - flast ) // fbar = f + theta * ( f - flast )
// flast_ = f_ // flast_ = f_
for ( int i = 0; i < n_; i++ ) for ( int i = 0; i < n_; i++ )
...@@ -64,6 +64,6 @@ void AndersonMixer::update(const double* f, double* theta, double* fbar) ...@@ -64,6 +64,6 @@ void AndersonMixer::update(const double* f, double* theta, double* fbar)
const double ftmp = f[i]; const double ftmp = f[i];
fbar[i] = ftmp + *theta * ( ftmp - flast_[i] ); fbar[i] = ftmp + *theta * ( ftmp - flast_[i] );
flast_[i] = ftmp; flast_[i] = ftmp;
} }
extrapolate_ = true; extrapolate_ = true;
} }
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// AndersonMixer.h // AndersonMixer.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: AndersonMixer.h,v 1.4 2007-03-17 01:14:00 fgygi Exp $ // $Id: AndersonMixer.h,v 1.5 2007-10-19 16:24:03 fgygi Exp $
#ifndef ANDERSONMIXER_H #ifndef ANDERSONMIXER_H
#define ANDERSONMIXER_H #define ANDERSONMIXER_H
...@@ -18,7 +18,7 @@ class AndersonMixer ...@@ -18,7 +18,7 @@ class AndersonMixer
const Context* const pctxt_; // pointer to relevant Context, null if local const Context* const pctxt_; // pointer to relevant Context, null if local
double theta_max_; // maximum extrapolation double theta_max_; // maximum extrapolation
double theta_nc_; // negative curvature value double theta_nc_; // negative curvature value
std::valarray<double> flast_; // last residual std::valarray<double> flast_; // last residual
bool extrapolate_; // state variable bool extrapolate_; // state variable
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// AngleCmd.h: // AngleCmd.h:
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: AngleCmd.h,v 1.1 2005-06-27 22:35:26 fgygi Exp $ // $Id: AngleCmd.h,v 1.2 2007-10-19 16:24:03 fgygi Exp $
#ifndef ANGLECMD_H #ifndef ANGLECMD_H
#define ANGLECMD_H #define ANGLECMD_H
...@@ -24,7 +24,7 @@ class AngleCmd : public Cmd ...@@ -24,7 +24,7 @@ class AngleCmd : public Cmd
char *name(void) const { return "angle"; } char *name(void) const { return "angle"; }
char *help_msg(void) const char *help_msg(void) const
{ {
return return
"\n angle\n\n" "\n angle\n\n"
" syntax: angle name1 name2 name3\n\n" " syntax: angle name1 name2 name3\n\n"
" The angle command prints the angle defined by three atoms.\n\n"; " The angle command prints the angle defined by three atoms.\n\n";
...@@ -40,15 +40,15 @@ class AngleCmd : public Cmd ...@@ -40,15 +40,15 @@ class AngleCmd : public Cmd
} }
return 1; return 1;
} }
string name1(argv[1]); string name1(argv[1]);
string name2(argv[2]); string name2(argv[2]);
string name3(argv[3]); string name3(argv[3]);
Atom* a1 = s->atoms.findAtom(name1); Atom* a1 = s->atoms.findAtom(name1);
Atom* a2 = s->atoms.findAtom(name2); Atom* a2 = s->atoms.findAtom(name2);
Atom* a3 = s->atoms.findAtom(name3); Atom* a3 = s->atoms.findAtom(name3);
if ( a1 == 0 || a2 == 0 || a3 == 0 ) if ( a1 == 0 || a2 == 0 || a3 == 0 )
{ {
if ( ui->onpe0() ) if ( ui->onpe0() )
{ {
if ( a1 == 0 ) if ( a1 == 0 )
...@@ -63,28 +63,28 @@ class AngleCmd : public Cmd ...@@ -63,28 +63,28 @@ class AngleCmd : public Cmd
} }
return 1; return 1;
} }
if ( a1 == a2 || a2 == a3 || a3 == a1 ) if ( a1 == a2 || a2 == a3 || a3 == a1 )
{ {
if ( ui->onpe0() ) if ( ui->onpe0() )
{ {
cout << " <!-- AngleCmd: replicated atoms in " << name1 cout << " <!-- AngleCmd: replicated atoms in " << name1
<< " " << name2 << " " << name3 << " -->" << endl; << " " << name2 << " " << name3 << " -->" << endl;
} }
return 1; return 1;
} }
D3vector r12(a1->position()-a2->position()); D3vector r12(a1->position()-a2->position());
D3vector r32(a3->position()-a2->position()); D3vector r32(a3->position()-a2->position());
if ( norm(r12) == 0.0 || norm(r32) == 0.0 ) if ( norm(r12) == 0.0 || norm(r32) == 0.0 )
{ {
if ( ui->onpe0() ) if ( ui->onpe0() )
{ {
cout << " <!-- AngleCmd: atoms are too close -->" << endl; cout << " <!-- AngleCmd: atoms are too close -->" << endl;
} }
return 1; return 1;
} }
const double sp = normalized(r12) * normalized(r32); const double sp = normalized(r12) * normalized(r32);
const double c = max(-1.0,min(1.0,sp)); const double c = max(-1.0,min(1.0,sp));
const double a = (180.0/M_PI)*acos(c); const double a = (180.0/M_PI)*acos(c);
...@@ -92,9 +92,9 @@ class AngleCmd : public Cmd ...@@ -92,9 +92,9 @@ class AngleCmd : public Cmd
{ {
cout.setf(ios::fixed,ios::floatfield); cout.setf(ios::fixed,ios::floatfield);
cout << " <!-- angle " << name1 << "-" << name2 << "-" << name3 cout << " <!-- angle " << name1 << "-" << name2 << "-" << name3
<< ": " << ": "
<< setprecision(3) << setprecision(3)
<< a << " (deg) -->" << endl; << a << " (deg) -->" << endl;
} }
return 0; return 0;
} }
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
// AngleConstraint.C // AngleConstraint.C
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: AngleConstraint.C,v 1.2 2005-09-16 23:08:11 fgygi Exp $ // $Id: AngleConstraint.C,v 1.3 2007-10-19 16:24:03 fgygi Exp $
#include "AngleConstraint.h" #include "AngleConstraint.h"
#include "AtomSet.h" #include "AtomSet.h"
...@@ -26,7 +26,7 @@ void AngleConstraint::setup(const AtomSet& atoms) ...@@ -26,7 +26,7 @@ void AngleConstraint::setup(const AtomSet& atoms)
m1_ = atoms.species_list[is1_]->mass() * 1822.89; m1_ = atoms.species_list[is1_]->mass() * 1822.89;
assert(m1_>0.0); assert(m1_>0.0);
m1_inv_ = 1.0 / m1_; m1_inv_ = 1.0 / m1_;
is2_ = atoms.is(name2_); is2_ = atoms.is(name2_);
ia2_ = atoms.ia(name2_); ia2_ = atoms.ia(name2_);
assert(is2_>=0); assert(is2_>=0);
...@@ -34,7 +34,7 @@ void AngleConstraint::setup(const AtomSet& atoms) ...@@ -34,7 +34,7 @@ void AngleConstraint::setup(const AtomSet& atoms)
m2_ = atoms.species_list[is2_]->mass() * 1822.89; m2_ = atoms.species_list[is2_]->mass() * 1822.89;
assert(m2_>0.0); assert(m2_>0.0);
m2_inv_ = 1.0 / m2_; m2_inv_ = 1.0 / m2_;
is3_ = atoms.is(name3_); is3_ = atoms.is(name3_);
ia3_ = atoms.ia(name3_); ia3_ = atoms.ia(name3_);
assert(is3_>=0); assert(is3_>=0);
...@@ -63,7 +63,7 @@ vector<vector<double> > &rp) const ...@@ -63,7 +63,7 @@ vector<vector<double> > &rp) const
double* pr1p = &rp[is1_][3*ia1_]; double* pr1p = &rp[is1_][3*ia1_];
double* pr2p = &rp[is2_][3*ia2_]; double* pr2p = &rp[is2_][3*ia2_];
double* pr3p = &rp[is3_][3*ia3_]; double* pr3p = &rp[is3_][3*ia3_];
D3vector r1(pr1); D3vector r1(pr1);
D3vector r2(pr2); D3vector r2(pr2);
D3vector r3(pr3); D3vector r3(pr3);
...@@ -72,7 +72,7 @@ vector<vector<double> > &rp) const ...@@ -72,7 +72,7 @@ vector<vector<double> > &rp) const
const double a = bond_angle(r1,r2,r3); const double a = bond_angle(r1,r2,r3);
double ng = g1*g1 + g2*g2 + g3*g3; double ng = g1*g1 + g2*g2 + g3*g3;
assert(ng>=0.0); assert(ng>=0.0);
D3vector r1p(pr1p); D3vector r1p(pr1p);
D3vector r2p(pr2p); D3vector r2p(pr2p);
D3vector r3p(pr3p); D3vector r3p(pr3p);
...@@ -81,7 +81,7 @@ vector<vector<double> > &rp) const ...@@ -81,7 +81,7 @@ vector<vector<double> > &rp) const
const double ap = bond_angle(r1p,r2p,r3p); const double ap = bond_angle(r1p,r2p,r3p);
const double ngp = g1p*g1p + g2p*g2p + g3p*g3p; const double ngp = g1p*g1p + g2p*g2p + g3p*g3p;
assert(ngp>=0.0); assert(ngp>=0.0);
const double err = fabs(ap-angle_); const double err = fabs(ap-angle_);
if ( ng == 0.0 ) if ( ng == 0.0 )
...@@ -119,25 +119,25 @@ vector<vector<double> > &rp) const ...@@ -119,25 +119,25 @@ vector<vector<double> > &rp) const
// no action necessary // no action necessary
} }
} }
// test alignment of the gradient at r and at rp // test alignment of the gradient at r and at rp
// compute scalar product of normalized gradients // compute scalar product of normalized gradients
const double sp = ( g1*g1p + g2*g2p + g3*g3p ) / sqrt( ng * ngp ); const double sp = ( g1*g1p + g2*g2p + g3*g3p ) / sqrt( ng * ngp );
if ( fabs(sp) < 0.5*sqrt(3.0) ) if ( fabs(sp) < 0.5*sqrt(3.0) )
{ {
g1 = g1p; g1 = g1p;
g2 = g2p; g2 = g2p;
g3 = g3p; g3 = g3p;
} }
const double den = m1_inv_ * g1 * g1p + const double den = m1_inv_ * g1 * g1p +
m2_inv_ * g2 * g2p + m2_inv_ * g2 * g2p +
m3_inv_ * g3 * g3p; m3_inv_ * g3 * g3p;
#if DEBUG_CONSTRAINTS #if DEBUG_CONSTRAINTS
cout << " <!-- "; cout << " <!-- ";
cout << " AngleConstraint::enforce_r: " cout << " AngleConstraint::enforce_r: "
<< name1_ << " " << name2_ << " " << name3_ << name1_ << " " << name2_ << " " << name3_
<< " angle = " << ap << endl; << " angle = " << ap << endl;
cout << " AngleConstraint::enforce_r: r1 = " << r1 << endl; cout << " AngleConstraint::enforce_r: r1 = " << r1 << endl;
...@@ -159,25 +159,25 @@ vector<vector<double> > &rp) const ...@@ -159,25 +159,25 @@ vector<vector<double> > &rp) const
cout << " -->" << endl; cout << " -->" << endl;
#endif #endif
if ( err < tol_ ) return true; if ( err < tol_ ) return true;
// make one shake iteration // make one shake iteration
assert(fabs(den)>0.0); assert(fabs(den)>0.0);
const double lambda = - sigma(r1p,r2p,r3p) / den; const double lambda = - sigma(r1p,r2p,r3p) / den;
pr1p[0] += m1_inv_ * lambda * g1.x; pr1p[0] += m1_inv_ * lambda * g1.x;
pr1p[1] += m1_inv_ * lambda * g1.y; pr1p[1] += m1_inv_ * lambda * g1.y;
pr1p[2] += m1_inv_ * lambda * g1.z; pr1p[2] += m1_inv_ * lambda * g1.z;
pr2p[0] += m2_inv_ * lambda * g2.x; pr2p[0] += m2_inv_ * lambda * g2.x;
pr2p[1] += m2_inv_ * lambda * g2.y; pr2p[1] += m2_inv_ * lambda * g2.y;
pr2p[2] += m2_inv_ * lambda * g2.z; pr2p[2] += m2_inv_ * lambda * g2.z;
pr3p[0] += m3_inv_ * lambda * g3.x; pr3p[0] += m3_inv_ * lambda * g3.x;
pr3p[1] += m3_inv_ * lambda * g3.y; pr3p[1] += m3_inv_ * lambda * g3.y;
pr3p[2] += m3_inv_ * lambda * g3.z; pr3p[2] += m3_inv_ * lambda * g3.z;
return false; return false;
} }
...@@ -191,30 +191,30 @@ vector<vector<double> > &v0) const ...@@ -191,30 +191,30 @@ vector<vector<double> > &v0) const
double* pv1 = &v0[is1_][3*ia1_]; double* pv1 = &v0[is1_][3*ia1_];
double* pv2 = &v0[is2_][3*ia2_]; double* pv2 = &v0[is2_][3*ia2_];
double* pv3 = &v0[is3_][3*ia3_]; double* pv3 = &v0[is3_][3*ia3_];
D3vector r1(pr1); D3vector r1(pr1);
D3vector r2(pr2); D3vector r2(pr2);
D3vector r3(pr3); D3vector r3(pr3);
D3vector g1,g2,g3; D3vector g1,g2,g3;
grad_sigma(r1,r2,r3,g1,g2,g3); grad_sigma(r1,r2,r3,g1,g2,g3);
D3vector v1(pv1); D3vector v1(pv1);
D3vector v2(pv2); D3vector v2(pv2);
D3vector v3(pv3); D3vector v3(pv3);
const double proj = v1*g1 + v2*g2 + v3*g3; const double proj = v1*g1 + v2*g2 + v3*g3;
const double norm2 = g1*g1 + g2*g2 + g3*g3; const double norm2 = g1*g1 + g2*g2 + g3*g3;
// if the gradient is too small, do not attempt correction // if the gradient is too small, do not attempt correction
if ( norm2 < 1.e-6 ) return true; if ( norm2 < 1.e-6 ) return true;
const double err = fabs(proj)/sqrt(norm2); const double err = fabs(proj)/sqrt(norm2);
#if DEBUG_CONSTRAINTS #if DEBUG_CONSTRAINTS
cout << " <!-- "; cout << " <!-- ";
cout << " AngleConstraint::enforce_v: " cout << " AngleConstraint::enforce_v: "
<< name1_ << " " << name2_ << " " << name3_ << endl; << name1_ << " " << name2_ << " " << name3_ << endl;
cout << " AngleConstraint::enforce_v: tol = " << tol_ << endl; cout << " AngleConstraint::enforce_v: tol = " << tol_ << endl;
cout << " AngleConstraint::enforce_v: err = " << err << endl; cout << " AngleConstraint::enforce_v: err = " << err << endl;
...@@ -224,27 +224,27 @@ vector<vector<double> > &v0) const ...@@ -224,27 +224,27 @@ vector<vector<double> > &v0) const
cout << " -->" << endl; cout << " -->" << endl;
#endif #endif
if ( err < tol_ ) return true; if ( err < tol_ ) return true;
// make one shake iteration // make one shake iteration
const double den = m1_inv_ * g1 * g1 + const double den = m1_inv_ * g1 * g1 +
m2_inv_ * g2 * g2 + m2_inv_ * g2 * g2 +
m3_inv_ * g3 * g3; m3_inv_ * g3 * g3;
assert(den>0.0); assert(den>0.0);
const double eta = -proj / den; const double eta = -proj / den;
pv1[0] += m1_inv_ * eta * g1.x; pv1[0] += m1_inv_ * eta * g1.x;
pv1[1] += m1_inv_ * eta * g1.y; pv1[1] += m1_inv_ * eta * g1.y;
pv1[2] += m1_inv_ * eta * g1.z; pv1[2] += m1_inv_ * eta * g1.z;
pv2[0] += m2_inv_ * eta * g2.x; pv2[0] += m2_inv_ * eta * g2.x;
pv2[1] += m2_inv_ * eta * g2.y; pv2[1] += m2_inv_ * eta * g2.y;
pv2[2] += m2_inv_ * eta * g2.z; pv2[2] += m2_inv_ * eta * g2.z;
pv3[0] += m3_inv_ * eta * g3.x; pv3[0] += m3_inv_ * eta * g3.x;
pv3[1] += m3_inv_ * eta * g3.y; pv3[1] += m3_inv_ * eta * g3.y;
pv3[2] += m3_inv_ * eta * g3.z; pv3[2] += m3_inv_ * eta * g3.z;
return false; return false;
} }
...@@ -258,30 +258,30 @@ void AngleConstraint::compute_force(const vector<vector<double> > &r0, ...@@ -258,30 +258,30 @@ void AngleConstraint::compute_force(const vector<vector<double> > &r0,
const double* pf1 = &f[is1_][3*ia1_]; const double* pf1 = &f[is1_][3*ia1_];
const double* pf2 = &f[is2_][3*ia2_]; const double* pf2 = &f[is2_][3*ia2_];
const double* pf3 = &f[is3_][3*ia3_]; const double* pf3 = &f[is3_][3*ia3_];
D3vector r1(pr1); D3vector r1(pr1);
D3vector r2(pr2); D3vector r2(pr2);
D3vector r3(pr3); D3vector r3(pr3);
D3vector g1,g2,g3; D3vector g1,g2,g3;
grad_sigma(r1,r2,r3,g1,g2,g3); grad_sigma(r1,r2,r3,g1,g2,g3);
D3vector f1(pf1); D3vector f1(pf1);
D3vector f2(pf2); D3vector f2(pf2);
D3vector f3(pf3); D3vector f3(pf3);
const double norm2 = g1*g1 + g2*g2 + g3*g3; const double norm2 = g1*g1 + g2*g2 + g3*g3;
if ( norm2 == 0.0 ) if ( norm2 == 0.0 )
{ {
force_ = 0.0; force_ = 0.0;
return; return;
} }
const double proj = f1*g1 + f2*g2 + f3*g3; const double proj = f1*g1 + f2*g2 + f3*g3;
force_ = -proj/norm2; force_ = -proj/norm2;
// compute weight // compute weight
const double z = m1_inv_ * g1 * g1 + const double z = m1_inv_ * g1 * g1 +
m2_inv_ * g2 * g2 + m2_inv_ * g2 * g2 +
m3_inv_ * g3 * g3; m3_inv_ * g3 * g3;
assert(z > 0.0); assert(z > 0.0);
...@@ -292,10 +292,10 @@ void AngleConstraint::compute_force(const vector<vector<double> > &r0, ...@@ -292,10 +292,10 @@ void AngleConstraint::compute_force(const vector<vector<double> > &r0,
const double r32s = norm(r3-r2); const double r32s = norm(r3-r2);
const double fac = 180.0/M_PI; const double fac = 180.0/M_PI;
const double cos_theta = normalized(r1-r2)*normalized(r3-r2); const double cos_theta = normalized(r1-r2)*normalized(r3-r2);
const double zcheck = fac*fac * const double zcheck = fac*fac *
( m1_inv_ / r12s + ( m1_inv_ / r12s +
m2_inv_ * ( (r12s+r32s-2*sqrt(r12s*r32s)*cos_theta) /(r12s*r32s) ) + m2_inv_ * ( (r12s+r32s-2*sqrt(r12s*r32s)*cos_theta) /(r12s*r32s) ) +
m3_inv_ / r32s m3_inv_ / r32s
); );
cout << " <!-- AngleConstraint: z=" << z << " zcheck=" << zcheck << " -->" cout << " <!-- AngleConstraint: z=" << z << " zcheck=" << zcheck << " -->"
<< endl; << endl;
...@@ -310,7 +310,7 @@ double AngleConstraint::sigma(D3vector a, D3vector b, D3vector c) const ...@@ -310,7 +310,7 @@ double AngleConstraint::sigma(D3vector a, D3vector b, D3vector c) const
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2, void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
const D3vector &r3, const D3vector &r3,
D3vector &g1, D3vector &g2, D3vector &g3) const D3vector &g1, D3vector &g2, D3vector &g3) const
{ {
...@@ -321,7 +321,7 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2, ...@@ -321,7 +321,7 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
D3vector e12(normalized(r12)); D3vector e12(normalized(r12));
D3vector e32(normalized(r32)); D3vector e32(normalized(r32));
const double ss = length(e12^e32); const double ss = length(e12^e32);
// use simplified expression if the angle is smaller than 0.2 degrees // use simplified expression if the angle is smaller than 0.2 degrees
if ( ss < sin(0.2*(M_PI/180.0) ) ) if ( ss < sin(0.2*(M_PI/180.0) ) )
{ {
...@@ -351,43 +351,43 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2, ...@@ -351,43 +351,43 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
{ {
// angle is large enough. Use finite differences // angle is large enough. Use finite differences
//cout << " ========= grad_sigma using finite differences" << endl; //cout << " ========= grad_sigma using finite differences" << endl;
const double r12_inv = 1.0/length(r12); const double r12_inv = 1.0/length(r12);
const double r32_inv = 1.0/length(r32); const double r32_inv = 1.0/length(r32);
const double a = bond_angle(r1,r2,r3); const double a = bond_angle(r1,r2,r3);
const double l12 = length(r1-r2); const double l12 = length(r1-r2);
const double l32 = length(r3-r2); const double l32 = length(r3-r2);
// displacement h causes changes in angle of less than 0.05 degrees // displacement h causes changes in angle of less than 0.05 degrees
const double h = 0.05 * (M_PI/180.0) * min(l12,l32); const double h = 0.05 * (M_PI/180.0) * min(l12,l32);
const double fac = 0.5 / h; const double fac = 0.5 / h;
D3vector dx(h,0,0), dy(0,h,0), dz(0,0,h); D3vector dx(h,0,0), dy(0,h,0), dz(0,0,h);
// compute gradient at r // compute gradient at r
g1.x = fac * ( sigma(r1+dx,r2,r3) - sigma(r1-dx,r2,r3) ); g1.x = fac * ( sigma(r1+dx,r2,r3) - sigma(r1-dx,r2,r3) );
g1.y = fac * ( sigma(r1+dy,r2,r3) - sigma(r1-dy,r2,r3) ); g1.y = fac * ( sigma(r1+dy,r2,r3) - sigma(r1-dy,r2,r3) );
g1.z = fac * ( sigma(r1+dz,r2,r3) - sigma(r1-dz,r2,r3) ); g1.z = fac * ( sigma(r1+dz,