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
22 #include "Space_warper.h"
23
24 #include "MatrixAlgebra.h"
25 #include "Sample_point.h"
26 #include "qmc_io.h"
27
Space_warper()28 Space_warper::Space_warper() {
29 weight_basis=NULL;
30 vector <string> basistxt;
31 basistxt.push_back("NULL");
32 basistxt.push_back("POLYPADE");
33 basistxt.push_back("RCUT");
34 basistxt.push_back("2.0");
35 basistxt.push_back("BETA0");
36 basistxt.push_back("-.9");
37 basistxt.push_back("NFUNC");
38 basistxt.push_back("1");
39 allocate(basistxt, weight_basis);
40 warp_on=1;
41 ex=1;
42 }
43
read(vector<string> & words)44 void Space_warper::read(vector <string> & words) {
45 unsigned int pos=0;
46 readvalue(words, pos=0, ex, "EXPONENT");
47 }
48
49
50 //----------------------------------------------------------------------
transformation_matrix(const Array2<doublevar> & ellePrime,const Array2<doublevar> & elle,Array2<doublevar> & ratio)51 void transformation_matrix (const Array2 <doublevar> &ellePrime,
52 const Array2 <doublevar> &elle,
53 Array2 <doublevar> &ratio) {
54
55 assert(elle.GetDim(0)==elle.GetDim(1)==3);
56 assert(ellePrime.GetDim(0)==ellePrime.GetDim(1)==3);
57 ratio.Resize(3,3);
58 Array2 <doublevar> elleInverse(3,3,0.0);
59
60 InvertMatrix(elle, elleInverse, 3);
61 MultiplyMatrices(ellePrime, elleInverse, ratio, 3);
62
63 }
64
65 //----------------------------------------------------------------------
space_warp(Sample_point * refsample,Sample_point * sample,int e,Array1<doublevar> & R_old,Array1<doublevar> & R_new,doublevar & jacobian)66 int Space_warper::space_warp (Sample_point * refsample, Sample_point * sample,
67 int e, Array1 <doublevar> & R_old,
68 Array1 <doublevar> & R_new,
69 doublevar & jacobian ) {
70
71 int natoms=refsample->ionSize();
72 int moved=0;
73 refsample->updateEIDist();
74 Array1 <doublevar> dist(5);
75
76 Array2 <doublevar> jacob_matrix(3,3);
77 jacob_matrix=0;
78 assert(R_old.GetDim(0)==3);
79 R_new.Resize(3);
80 R_new=R_old;
81
82
83 if(!warp_on)
84 { jacobian=1; return 0; }
85
86 //Here we stretch or shrink the cell appropriately to match the
87 //different lattice parameters. This is a guaranteed way to
88 //make sure the two spaces actually have a 1-1 relationship in a PBC calculation.
89 Array2 <doublevar> reflatvec;
90 Array1 <doublevar> reforigin;
91 if(refsample->getBounds(reflatvec, reforigin)) {
92 Array2 <doublevar> latvec;
93 Array1 <doublevar> origin;
94 sample->getBounds(latvec, origin);
95 assert(reflatvec.GetDim(0)==3 && reflatvec.GetDim(1)==3 && reforigin.GetDim(0)==3);
96 assert(latvec.GetDim(0)==3 && latvec.GetDim(1)==3 && origin.GetDim(0)==3);
97
98 transformation_matrix(latvec, reflatvec, jacob_matrix);
99 R_new=0;
100 //cout << " jacobian " << endl;
101 for(int i=0; i< 3; i++) {
102 for(int j=0; j< 3; j++) {
103 R_new(i)+=jacob_matrix(i,j)*R_old(j);
104 //cout << jacob_matrix(i,j) << " ";
105 }
106 //cout << endl;
107 }
108
109 //cout << "rnew rold \n";
110 //for(int i=0; i< 3; i++) cout << R_new(i) <<" " << R_old(i) << endl;
111 //cout << endl;
112 jacobian=Determinant(jacob_matrix,3);
113 return 1;
114 }
115
116
117 //If it is not a periodic calculation, then space-warp using
118 //either a general basis function or the 1/r warping suggested
119 //by Filippi and Umrigar(using 1/r by default because the full
120 //estimator works better with it)
121
122 //doublevar cutoff=weight_basis->cutoff(0);
123 doublevar cutoff=1e99;
124
125 //Array2 <doublevar> basisval(1,5);
126
127 Array1 <doublevar> displace(3, 0.0);
128 //if(0) {
129 doublevar norm=0;
130 Array1 <doublevar> der_norm(3,0.0);
131 Array1 <doublevar> weights(natoms,0.0);
132 Array2 <doublevar> weight_der(natoms, 3,0.0);
133 Array2 <doublevar> deltar(natoms, 3,0.0);
134
135 Array1 <doublevar> oldrefpos(3);
136 refsample->getElectronPos(e,oldrefpos);
137
138
139
140
141 if(natoms >0) {
142 refsample->setElectronPosNoNotify(e,R_old);
143 refsample->updateEIDist();
144 }
145
146 for(int at=0; at < natoms; at++) {
147
148 refsample->getEIDist(e,at, dist);
149 //cout << "atom " << at << "dist " << dist(0) << endl;
150 if(dist(0) < cutoff) {
151 //if(moved==1) cout << "Warning--warped twice " << endl;
152
153 moved=1;
154 Array1 <doublevar> refionpos(3);
155 refsample->getIonPos(at, refionpos);
156 Array1 <doublevar> newionpos(3);
157 sample->getIonPos(at, newionpos);
158
159 //weight_basis->calcLap(dist, basisval);
160 doublevar z=1.0;//refsample->getIonCharge(at);
161
162
163
164 //doublevar func=z/(dist(1)*dist(1));
165 //doublevar func=z/(dist(1));
166 doublevar func=z/(pow(dist(0),ex));
167 //doublevar func=basisval(0,0);
168 norm+=func;
169
170 for(int d=0; d< 3; d++) {
171 displace(d) += func*(newionpos(d)-refionpos(d));
172 }
173
174 weights(at)=func;
175 for(int d=0; d< 3; d++) {
176 deltar(at,d)=newionpos(d)-refionpos(d);
177 //weight_der(at,d)= -4*z*dist(d+2)/(dist(1)*dist(1)*dist(1));
178 weight_der(at,d)= -ex*z*dist(d+2)/pow(dist(0),ex+2);
179 //weight_der(at,d)=basisval(0,d+1);
180 }
181 for(int d=0; d< 3; d++) der_norm(d)+=weight_der(at,d);
182 }
183 }
184
185
186 if(natoms >0) {
187 refsample->setElectronPosNoNotify(e,oldrefpos);
188 refsample->updateEIDist();
189 }
190
191 //cout << "norm " << norm << endl;
192 if(moved && norm > 1e-14) {
193 for(int d=0; d< 3; d++) {
194 R_new(d)+=displace(d)/norm;
195 }
196
197 jacob_matrix=0;
198 for(int i=0; i< 3; i++)
199 for(int j=0; j< 3; j++)
200 for(int at=0; at < natoms; at++) {
201 jacob_matrix(i,j)+=deltar(at,i)*(weight_der(at,j)/norm
202 -weights(at)*der_norm(j)/(norm*norm));
203 //cout << " at " << at << "der " << weight_der(at,j) << " weights "
204 // << weights(at) << " der_norm " << der_norm(j)
205 // << endl;
206 }
207 }
208 //}
209
210
211 for(int i=0; i< 3; i++) jacob_matrix(i,i)+=1;
212
213 //cout << "jacobian matrix " << endl;
214 //for(int i=0; i< 3; i++) {
215 // for(int j=0; j< 3; j++)
216 // cout << jacob_matrix(i,j) << " ";
217 // cout << endl;
218 //}
219 jacobian=Determinant(jacob_matrix, 3);
220 //cout << "jacobian " << jacobian << endl;
221 return moved;
222 }
223
224 //----------------------------------------------------------------------
225
warp_all(Sample_point * refsample,Sample_point * sample)226 doublevar Space_warper::warp_all(Sample_point * refsample,Sample_point * sample) {
227 doublevar jacobian=1,tot_jacobian=1;
228 int nelectrons=refsample->electronSize();
229 for(int e=0; e< nelectrons; e++) {
230 Array1 <doublevar> ref_pos(3);
231 Array1 <doublevar> warp_pos(3);
232 refsample->getElectronPos(e,ref_pos);
233 space_warp(refsample,sample, e, ref_pos, warp_pos, jacobian);
234 tot_jacobian*=jacobian;
235 sample->setElectronPos(e,warp_pos);
236 }
237 return tot_jacobian;
238 }
239