1 /*
2
3 Copyright (C) 2007 Lucas K. Wagner
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18
19 */
20
21 #include "Force_fitter.h"
22 #include "MatrixAlgebra.h"
23
setup(doublevar maxcut,int nexpansion)24 void Force_fitter::setup(doublevar maxcut, int nexpansion ) {
25 m=2;
26 cut=maxcut;
27 nexp=nexpansion;
28 coeff.Resize(nexp);
29 Array2 <doublevar> S(nexp, nexp);
30 Array1 <doublevar> h(nexp);
31 for(int i=0; i< nexp; i++) {
32 for(int j=0; j< nexp; j++) {
33 S(i,j)=pow(cut,m+i+j+1)/(m+i+j+1);
34 }
35 h(i)=pow(cut,i+1)/(i+1);
36 }
37
38 Array2 <doublevar> Sinv(nexp, nexp);
39 InvertMatrix(S, Sinv, nexp);
40
41 coeff=0;
42 for(int i=0; i< nexp; i++)
43 for(int j=0;j< nexp;j++)
44 coeff(i)+=Sinv(i,j)*h(j);
45 }
46
47
fit_force(Array1<doublevar> & r,Array1<doublevar> & bare_force,Array1<doublevar> & fit_force)48 void Force_fitter::fit_force(Array1 <doublevar> & r,
49 Array1 <doublevar> & bare_force,
50 Array1 <doublevar> & fit_force) {
51 int ndim=bare_force.GetDim(0);
52
53 fit_force=bare_force; return;
54
55 fit_force.Resize(ndim);
56 fit_force=0;
57
58 if(r(0) < cut) {
59 for(int d=0; d< ndim; d++) {
60 doublevar basis=0;
61 for(int i=0; i< nexp; i++) {
62 basis+=coeff(i)*pow(r(0),i+m);
63 }
64 fit_force(d)=bare_force(d)*basis;
65 }
66 }
67 else {
68 fit_force=bare_force;
69 }
70
71
72 }
73
74