cutoff.cpp 2.59 KB
Newer Older
1
#include "soap/cutoff.hpp"
Carl Poelking's avatar
Carl Poelking committed
2
3
4
5
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
Carl Poelking's avatar
Carl Poelking committed
6
7
8
9
10
11
12
13

namespace soap {


void CutoffFunction::configure(Options &options) {
    _Rc = options.get<double>("radialcutoff.Rc");
    _Rc_width = options.get<double>("radialcutoff.Rc_width");
    _center_weight = options.get<double>("radialcutoff.center_weight");
Carl Poelking's avatar
Carl Poelking committed
14
15
16
17

    GLOG() << "Weighting function with "
    	<< "Rc = " << _Rc
    	<< ", _Rc_width = " << _Rc_width
18
		<< ", central weight = " << _center_weight << std::endl;
Carl Poelking's avatar
Carl Poelking committed
19
20
}

Carl Poelking's avatar
Carl Poelking committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
bool CutoffFunction::isWithinCutoff(double r) {
    return (r <= _Rc);
}

double CutoffFunction::calculateWeight(double r) {
    double weight_at_r = 1.;
    if (r > _Rc) {
        weight_at_r = -1.e-10;
    }
    else if (r <= _Rc - _Rc_width) {
        weight_at_r = 1.;
    }
    else {
        weight_at_r = 0.5*(1+cos(M_PI*(r-_Rc+_Rc_width)/_Rc_width));
        //weight_at_r = pow(cos(0.5*M_PI*(r-_Rc+_Rc_width)/_Rc_width),2);
    }
    return weight_at_r;
}

vec CutoffFunction::calculateGradientWeight(double r, vec d) {
    vec grad_weight(0.,0.,0.);
42
    if (r > _Rc || r <= _Rc - _Rc_width) {
Carl Poelking's avatar
Carl Poelking committed
43
44
45
46
47
48
49
50
51
52
        ;
    }
    else {
        grad_weight = -0.5*sin(M_PI*(r-_Rc+_Rc_width)/_Rc_width)*M_PI/_Rc_width * d;
        //double phi = 0.5*M_PI*(r-_Rc+_Rc_width)/_Rc_width;
        //grad_weight = - M_PI/_Rc_width*sin(phi)*cos(phi) * d;
    }
    return grad_weight;
}

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
void CutoffFunctionHeaviside::configure(Options &options) {
    _Rc = options.get<double>("radialcutoff.Rc");
    _Rc_width = options.get<double>("radialcutoff.Rc_width");
    _center_weight = options.get<double>("radialcutoff.center_weight");

    // To harmonize with adaptive Gaussian basis sets, increase cutoff:
    if (options.hasKey("radialcutoff.Rc_heaviside")) {
        _Rc = options.get<double>("radialcutoff.Rc_heaviside"); // <- Set by RadialBasisGaussian::configure
    }

    GLOG() << "Weighting function with "
        << "Rc = " << _Rc << ", central weight = " << _center_weight << std::endl;
}

bool CutoffFunctionHeaviside::isWithinCutoff(double r) {
    return (r <= _Rc);
}

double CutoffFunctionHeaviside::calculateWeight(double r) {
    return (r > _Rc) ? -1.e-10 : 1.;
}

vec CutoffFunctionHeaviside::calculateGradientWeight(double r, vec d) {
    return vec(0.,0.,0.);
}

Carl Poelking's avatar
Carl Poelking committed
79
80
void CutoffFunctionFactory::registerAll(void) {
	CutoffFunctionOutlet().Register<CutoffFunction>("shifted-cosine");
81
	CutoffFunctionOutlet().Register<CutoffFunctionHeaviside>("heaviside");
Carl Poelking's avatar
Carl Poelking committed
82
83
84
85
}

}

Carl Poelking's avatar
Carl Poelking committed
86
87
BOOST_CLASS_EXPORT_IMPLEMENT(soap::CutoffFunctionHeaviside);