1 /*
2
3 Copyright (C) 2007 Lucas K. Wagner, 2008 Jindrich Kolorenc
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 "Qmc_std.h"
22 #include "BCS_wf.h"
23 #include "MatrixAlgebra.h"
24 #include "Sample_point.h"
25 #include "BCS_wf_data.h"
26
27 //----------------------------------------------------------------------
28
~BCS_wf()29 BCS_wf::~BCS_wf() {
30 if(jastcof) delete jastcof;
31 }
32
generateStorage(Wavefunction_storage * & wfstore)33 void BCS_wf::generateStorage(Wavefunction_storage * & wfstore)
34 {
35 wfstore=new BCS_wf_storage;
36 BCS_wf_storage * store;
37 recast(wfstore, store);
38 jast.generateStorage(store->jast_store);
39 store->inverse=inverse;
40 store->detVal=detVal;
41 //store->moVal.Resize(moVal.GetDim(1),moVal.GetDim(2));
42 store->derivatives.Resize(derivatives.GetDim(0),derivatives.GetDim(2));
43 store->derivatives_2.Resize(derivatives.GetDim(0),derivatives.GetDim(2));
44 }
45
46
47 //----------------------------------------------------------------------
48
init(Wavefunction_data * wfdata)49 void BCS_wf::init(Wavefunction_data * wfdata)
50 {
51
52 BCS_wf_data * dataptr;
53 recast(wfdata, dataptr);
54 recast(wfdata, parent);
55
56 //nmo=dataptr->molecorb->getNmo();
57 ndet=1;
58 nelectrons.Resize(2);
59 nelectrons=dataptr->nelectrons;
60
61 int tote=nelectrons(0)+nelectrons(1);
62 ndim=3;
63
64 if(nelectrons(0) != nelectrons(1)) {
65 error("For BCS wave function, the number of spin up electrons must"
66 " be equal to the number of spin down electrons.");
67 }
68
69 spin.Resize(tote);
70 spin=dataptr->spin;
71
72 //moVal.Resize(tote, nmo,5);
73 //updatedMoVal.Resize(nmo,5);
74
75 detVal.Resize (ndet);
76 inverse.Resize(ndet);
77
78 for(int det=0; det < ndet; det++) {
79 inverse(det).Resize(nelectrons(0), nelectrons(0));
80 inverse(det)=0;
81 for(int e=0; e< nelectrons(0); e++) {
82 inverse(det)(e,e)=1;
83 inverse(det)(e,e)=1;
84 }
85
86 detVal(det)=1;
87 }
88
89
90 electronIsStaleVal.Resize(tote);
91 electronIsStaleLap.Resize(tote);
92
93 electronIsStaleVal=0;
94 electronIsStaleLap=0;
95 updateEverythingVal=1;
96 updateEverythingLap=1;
97 sampleAttached=0;
98 dataAttached=0;
99
100 jast.init(&parent->jastdata);
101
102 derivatives.Resize(2*nelectrons(0),nelectrons(1),4);
103
104 assert(jastcof==NULL);
105 jastcof=new BCS_jastrow_u;
106
107 }
108
109 //----------------------------------------------------------------------
110
111 /*!
112 Behavior under staticSample: if sample_static is set, then
113 we ignore all electron moves in the main algorithm, and the only way
114 to update based on a move is by the getParmDepVal function. The
115 regular update functions will not work in this case.
116 */
notify(change_type change,int num)117 void BCS_wf::notify(change_type change, int num)
118 {
119
120 jast.notify(change,num);
121 switch(change)
122 {
123 case electron_move:
124 electronIsStaleVal(num)=1;
125 electronIsStaleLap(num)=1;
126 break;
127 case all_electrons_move:
128 updateEverythingVal=1;
129 updateEverythingLap=1;
130 break;
131 case wf_parm_change:
132 case all_wf_parms_change:
133 updateEverythingVal=1;
134 updateEverythingLap=1;
135 break;
136 case sample_attach:
137 sampleAttached=1;
138 updateEverythingVal=1;
139 updateEverythingLap=1;
140 break;
141 case data_attach:
142 dataAttached=1;
143 updateEverythingVal=1;
144 updateEverythingLap=1;
145 break;
146 default:
147 updateEverythingVal=1;
148 updateEverythingLap=1;
149 }
150
151 }
152
153 //----------------------------------------------------------------------
154
155
saveUpdate(Sample_point * sample,int e,Wavefunction_storage * wfstore)156 void BCS_wf::saveUpdate(Sample_point * sample, int e,
157 Wavefunction_storage * wfstore)
158 {
159 BCS_wf_storage * store;
160 recast(wfstore, store);
161 jast.saveUpdate(sample,e,store->jast_store);
162 store->inverse=inverse;
163 store->detVal=detVal;
164
165 //
166 // unpaired electrons go into extra one-body orbitals, this
167 // functionality will be resurrected later
168 //
169 //int nmo=moVal.GetDim(1);
170 //int nd=moVal.GetDim(2);
171 //for(int m=0; m< nmo; m++) {
172 // for(int d=0; d< nd; d++) {
173 // store->moVal(m,d)=moVal(e,m,d);
174 // }
175 //}
176
177 if(spin(e)==0) {
178 int nj=derivatives.GetDim(1);
179 assert(nj==nelectrons(0));
180 int nd=derivatives.GetDim(2);
181 for(int j=0; j< nj; j++) {
182 for(int d=0; d< nd; d++) {
183 store->derivatives(j,d)=derivatives(e,j,d);
184 }
185 }
186 int ep=e+nelectrons(0);
187 for(int j=0; j< nj; j++) {
188 for(int d=0; d< nd; d++) {
189 store->derivatives(j+nelectrons(0),d)=derivatives(ep,j,d);
190 }
191 }
192 }
193 else {
194 int ni=derivatives.GetDim(0);
195 assert(ni==2*nelectrons(0));
196 int nd=derivatives.GetDim(2);
197 int rede=e-nelectrons(0);
198 for(int i=0; i< ni; i++) {
199 for(int d=0; d< nd; d++) {
200 store->derivatives(i,d)=derivatives(i,rede,d);
201 }
202 }
203 }
204
205 }
206
207 //----------------------------------------------------------------------
208
restoreUpdate(Sample_point * sample,int e,Wavefunction_storage * wfstore)209 void BCS_wf::restoreUpdate(Sample_point * sample, int e,
210 Wavefunction_storage * wfstore)
211 {
212
213 BCS_wf_storage * store;
214 recast(wfstore, store);
215
216 jast.restoreUpdate(sample,e,store->jast_store);
217 detVal=store->detVal;
218 inverse=store->inverse;
219
220 //
221 // unpaired electrons go into extra one-body orbitals, this
222 // functionality will be resurrected later
223 //
224 //int nmo=moVal.GetDim(1);
225 //int nd=moVal.GetDim(2);
226 //for(int m=0; m< nmo; m++) {
227 // for(int d=0; d< nd; d++) {
228 // moVal(e,m,d)=store->moVal(m,d);
229 // }
230 //}
231
232 if(spin(e)==0) {
233 int nj=derivatives.GetDim(1);
234 assert(nj==nelectrons(0));
235 int nd=derivatives.GetDim(2);
236 for(int j=0; j< nj; j++) {
237 for(int d=0; d< nd; d++) {
238 derivatives(e,j,d)=store->derivatives(j,d);
239 }
240 }
241 int ep=e+nelectrons(0);
242 for(int j=0; j< nj; j++) {
243 for(int d=0; d< nd; d++) {
244 derivatives(ep,j,d)=store->derivatives(j+nelectrons(0),d);
245 }
246 }
247 }
248 else {
249 int ni=derivatives.GetDim(0);
250 assert(ni==2*nelectrons(0));
251 int nd=derivatives.GetDim(2);
252 int rede=e-nelectrons(0);
253 for(int i=0; i< ni; i++) {
254 for(int d=0; d< nd; d++) {
255 derivatives(i,rede,d)=store->derivatives(i,d);
256 }
257 }
258 }
259
260
261 electronIsStaleLap=0;
262 electronIsStaleVal=0;
263 updateEverythingLap=0;
264 updateEverythingVal=0;
265
266 //calcLap(sample);
267 }
268
269 //----------------------------------------------------------------------
270 //Added by Matous
271
saveUpdate(Sample_point * sample,int e1,int e2,Wavefunction_storage * wfstore)272 void BCS_wf::saveUpdate(Sample_point * sample, int e1, int e2,
273 Wavefunction_storage * wfstore)
274 {
275 BCS_wf_storage * store;
276 recast(wfstore, store);
277 jast.saveUpdate(sample,e1,e2,store->jast_store);
278 store->inverse=inverse;
279 store->detVal=detVal;
280
281 //
282 // unpaired electrons go into extra one-body orbitals, this
283 // functionality will be resurrected later
284 //
285 //int nmo=moVal.GetDim(1);
286 //int nd=moVal.GetDim(2);
287 //for(int m=0; m< nmo; m++) {
288 // for(int d=0; d< nd; d++) {
289 // store->moVal(m,d)=moVal(e,m,d);
290 // }
291 //}
292
293 // We assume the two electrons have opposite spins
294 assert( spin(e1) != spin(e2) );
295
296 int e_up, e_down;
297
298 e_up = spin(e1) ? e2 : e1;
299 e_down = spin(e1) ? e1 : e2;
300
301
302 {
303 int nj=derivatives.GetDim(1);
304 assert(nj==nelectrons(0));
305 int nd=derivatives.GetDim(2);
306 for(int j=0; j< nj; j++) {
307 for(int d=0; d< nd; d++) {
308 store->derivatives(j,d)=derivatives(e_up,j,d);
309 }
310 }
311 int ep=e_up+nelectrons(0);
312 for(int j=0; j< nj; j++) {
313 for(int d=0; d< nd; d++) {
314 store->derivatives(j+nelectrons(0),d)=derivatives(ep,j,d);
315 }
316 }
317 }
318 {
319 int ni=derivatives.GetDim(0);
320 assert(ni==2*nelectrons(0));
321 int nd=derivatives.GetDim(2);
322 int rede=e_down-nelectrons(0);
323 for(int i=0; i< ni; i++) {
324 for(int d=0; d< nd; d++) {
325 store->derivatives_2(i,d)=derivatives(i,rede,d);
326 }
327 }
328 }
329
330 }
331
332 //----------------------------------------------------------------------
333
restoreUpdate(Sample_point * sample,int e1,int e2,Wavefunction_storage * wfstore)334 void BCS_wf::restoreUpdate(Sample_point * sample, int e1, int e2,
335 Wavefunction_storage * wfstore)
336 {
337
338 BCS_wf_storage * store;
339 recast(wfstore, store);
340
341 jast.restoreUpdate(sample,e1,e2,store->jast_store);
342
343 detVal=store->detVal;
344 inverse=store->inverse;
345
346 //
347 // unpaired electrons go into extra one-body orbitals, this
348 // functionality will be resurrected later
349 //
350 //int nmo=moVal.GetDim(1);
351 //int nd=moVal.GetDim(2);
352 //for(int m=0; m< nmo; m++) {
353 // for(int d=0; d< nd; d++) {
354 // moVal(e,m,d)=store->moVal(m,d);
355 // }
356 //}
357
358 // We assume the two electrons have opposite spins
359 assert( spin(e1) != spin(e2) );
360
361 int e_up, e_down;
362
363 e_up = spin(e1) ? e2 : e1;
364 e_down = spin(e1) ? e1 : e2;
365
366 {
367 int ni=derivatives.GetDim(0);
368 assert(ni==2*nelectrons(0));
369 int nd=derivatives.GetDim(2);
370 int rede=e_down-nelectrons(0);
371 for(int i=0; i< ni; i++) {
372 for(int d=0; d< nd; d++) {
373 derivatives(i,rede,d)=store->derivatives_2(i,d);
374 }
375 }
376 }
377 {
378 int nj=derivatives.GetDim(1);
379 assert(nj==nelectrons(0));
380 int nd=derivatives.GetDim(2);
381 for(int j=0; j< nj; j++) {
382 for(int d=0; d< nd; d++) {
383 derivatives(e_up,j,d)=store->derivatives(j,d);
384 }
385 }
386 int ep=e_up+nelectrons(0);
387 for(int j=0; j< nj; j++) {
388 for(int d=0; d< nd; d++) {
389 derivatives(ep,j,d)=store->derivatives(j+nelectrons(0),d);
390 }
391 }
392 }
393
394
395 electronIsStaleLap=0;
396 electronIsStaleVal=0;
397 updateEverythingLap=0;
398 updateEverythingVal=0;
399
400
401 //calcLap(sample);
402 }
403
404 //----------------------------------------------------------------------
405
updateVal(Wavefunction_data * wfdata,Sample_point * sample)406 void BCS_wf::updateVal(Wavefunction_data * wfdata,
407 Sample_point * sample)
408 {
409 updateLap(wfdata,sample);
410 return;
411
412 }
413
414 //----------------------------------------------------------------------
415
updateLap(Wavefunction_data * wfdata,Sample_point * sample)416 void BCS_wf::updateLap( Wavefunction_data * wfdata,
417 Sample_point * sample)
418 {
419 assert(sampleAttached);
420 assert(dataAttached);
421
422
423
424 if(updateEverythingLap==1) {
425 calcLap(sample);
426 updateEverythingVal=0;
427 updateEverythingLap=0;
428 electronIsStaleLap=0;
429 electronIsStaleVal=0;
430 }
431 else {
432 for(int e=0; e< nelectrons(0)+nelectrons(1); e++) {
433 if(electronIsStaleLap(e)) {
434 if ( abs(detVal(0))>0 ) {
435 updateLap(e,sample);
436 } else {
437 calcLap(sample);
438 }
439 electronIsStaleLap(e)=0;
440 electronIsStaleVal(e)=0;
441 }
442 }
443
444 }
445
446
447 }
448
449 //-----------------------------------------------------------------------
450
451
getVal(Wavefunction_data * wfdata,int e,Wf_return & val)452 void BCS_wf::getVal(Wavefunction_data * wfdata, int e,
453 Wf_return & val)
454 {
455 assert(val.amp.GetDim(0) >=1);
456 assert(val.amp.GetDim(1) >= 1);
457
458 doublevar funcval=0;
459
460 for(int det=0; det < ndet; det++) {
461 funcval += detVal(det);
462 }
463 if(fabs(funcval) > 0)
464 val.amp(0,0)=log(fabs(funcval));
465 else val.amp(0,0)=-1e3;
466 double si=sign(funcval);
467
468 if(si>0)
469 val.phase(0,0)=0;
470 else val.phase(0,0)=pi;
471
472 }
473
474 //-------------------------------------------------------------------------
getSymmetricVal(Wavefunction_data * wfdata,int e,Wf_return & val)475 void BCS_wf::getSymmetricVal(Wavefunction_data * wfdata,
476 int e, Wf_return & val){
477
478 Array1 <doublevar> si(1, 0.0);
479 Array2 <doublevar> vals(1,1,0.0);
480 val.setVals(vals, si);
481 }
482
483
484 //----------------------------------------------------------------------
485
486
487 //Note that the matrix is almost symmetric--
488 //we can use that to halve the computational time
calcLap(Sample_point * sample)489 void BCS_wf::calcLap(Sample_point * sample)
490 {
491
492 int tote=nelectrons(0)+nelectrons(1);
493 Array3 <doublevar> temp_der;
494 Array2 <doublevar> temp_lap;
495 for(int e=0; e< tote; e++) {
496 sample->updateEIDist();
497 //
498 //parent->molecorb->updateLap(sample,e,0,updatedMoVal);
499 //for(int i=0; i< updatedMoVal.GetDim(0); i++) {
500 // for(int d=0; d< 5; d++) {
501 // moVal(e,i,d)=updatedMoVal(i,d);
502 // }
503 //}
504 }
505
506 Array2 <doublevar> onebody;
507 jast.updateLap(&parent->jastdata,sample);
508 jast.get_twobody(twobody);
509 jast.get_onebody_save(onebody);
510
511 int maxmatsize=max(nelectrons(0),nelectrons(1));
512 Array2 <doublevar> modet(maxmatsize, maxmatsize);
513 modet=0;
514 for(int det=0; det < ndet; det++ ) {
515 for(int i=0; i< nelectrons(0); i++) {
516 for(int j=0; j< nelectrons(1); j++) {
517 int jshift=nelectrons(0)+j;
518
519 Array2 <doublevar> ugrad;
520 jastcof->valGradLap(twobody,onebody,i,jshift,ugrad);
521 modet(i,j)=parent->magnification_factor*ugrad(0,0);
522
523 for(int d=0; d< 4; d++) {
524 derivatives(i,j,d)=ugrad(0,d+1);
525 derivatives(i+nelectrons(0),j,d)=ugrad(1,d+1);
526 }
527
528 }
529 //
530 // unpaired electrons go into extra one-body orbitals, this
531 // functionality will be resurrected later
532 //
533 //for(int j=nelectrons(1); j< nelectrons(0); j++) {
534 // modet(i,j)=moVal(i,parent->occupation(det,0)(j),0);
535 //}
536 }
537
538
539 detVal(det)=
540 TransposeInverseMatrix(modet,inverse(det), nelectrons(0)).val();
541
542 }
543
544 }
545
546 //------------------------------------------------------------------------
547
updateLap(int e,Sample_point * sample)548 void BCS_wf::updateLap(int e, Sample_point * sample ) {
549
550 int s=spin(e);
551 sample->updateEIDist();
552 //
553 //parent->molecorb->updateLap(sample,e,0,updatedMoVal);
554 //for(int i=0; i< updatedMoVal.GetDim(0); i++) {
555 // for(int d=0; d< 5; d++) {
556 // moVal(e,i,d)=updatedMoVal(i,d);
557 // }
558 //}
559
560 Array2 <doublevar> onebody;
561 jast.updateLap(&parent->jastdata,sample);
562 jast.get_twobody(twobody);
563 jast.get_onebody_save(onebody);
564
565
566 for(int det=0; det < ndet; det++ ) {
567 Array1 <doublevar> detUpdate(nelectrons(0));
568
569 if(s==0) { //------------spin up electron
570 for(int j=0; j< nelectrons(1); j++) {
571 int jshift=nelectrons(0)+j;
572 Array2 <doublevar> ugrad;
573 jastcof->valGradLap(twobody,onebody,e,jshift,ugrad);
574 detUpdate(j)=parent->magnification_factor*ugrad(0,0);
575
576 for(int d=0; d< 4; d++) {
577 derivatives(e,j,d)=ugrad(0,d+1);
578 derivatives(e+nelectrons(0),j,d)=ugrad(1,d+1);
579 }
580 }
581 //
582 // unpaired electrons go into extra one-body orbitals, this
583 // functionality will be resurrected later
584 //
585 //for(int j=nelectrons(1); j< nelectrons(0); j++) {
586 // detUpdate(j)=moVal(e,parent->occupation(det,0)(j),0);
587 //}
588 detVal(det)/=InverseUpdateColumn(inverse(det),detUpdate,e,nelectrons(0));
589
590 }
591 else { //--------------------spin down electrons
592 int eshift=e;
593 int edown=e-nelectrons(0);
594
595 for(int i=0; i< nelectrons(0); i++) {
596 Array2 <doublevar> ugrad;
597 jastcof->valGradLap(twobody,onebody,i,eshift,ugrad);
598 detUpdate(i)=parent->magnification_factor*ugrad(0,0);
599
600 for(int d=0; d< 4; d++) {
601 derivatives(i,edown,d)=ugrad(0,d+1);
602 derivatives(i+nelectrons(0),edown,d)=ugrad(1,d+1);
603 }
604 }
605 detVal(det)/=InverseUpdateRow(inverse(det),detUpdate,
606 edown,nelectrons(0));
607
608 }
609
610 }
611
612 }
613
614 //----------------------------------------------------------------------
615
616 /*!
617 */
618
getLap(Wavefunction_data * wfdata,int e,Wf_return & lap)619 void BCS_wf::getLap(Wavefunction_data * wfdata,
620 int e, Wf_return & lap)
621 {
622
623 getVal(wfdata,e,lap);
624
625 for(int d=1; d< 5; d++) {
626 if(spin(e)==0) {
627 for(int det=0; det < ndet; det++) {
628 doublevar temp=0;
629 for(int j=0; j< nelectrons(1); j++) {
630 temp+=derivatives(e,j,d-1)
631 *inverse(det)(e,j);
632 }
633 //
634 // unpaired electrons go into extra one-body orbitals, this
635 // functionality will be resurrected later
636 //
637 //for(int j=nelectrons(1); j < nelectrons(0); j++) {
638 // //error("need to do spin-polarized derivatives");
639 // temp+=moVal(e,parent->occupation(det,0)(j),d)
640 // *inverse(det)(e,j);
641 //}
642 if(ndet > 1) error("update BCS::getLap() for ndet > 1");
643 lap.amp(0,d)=parent->magnification_factor*temp;
644 }
645 }
646 else {
647 for(int det=0; det < ndet; det++) {
648 doublevar temp=0;
649 int red=parent->rede(e);
650 for(int i=0; i< nelectrons(0); i++) {
651 temp+=derivatives(i+nelectrons(0),red,d-1)
652 *inverse(det)(i,red);
653 }
654 lap.amp(0,d)=parent->magnification_factor*temp;
655 }
656
657 }
658 }
659
660 for(int d=1; d< 5; d++) {
661 lap.phase(0,d)=0.0;
662 }
663
664 }
665
666 //-------------------------------------------------------------------------
667
getParmDeriv(Wavefunction_data * wfdata,Sample_point * sample,Parm_deriv_return & parm_derivatives)668 int BCS_wf::getParmDeriv(Wavefunction_data *wfdata , Sample_point * sample,
669 Parm_deriv_return & parm_derivatives ) {
670
671 // analytic expressions below are valid only if jastrow can provide
672 // analytic derivatives (not sure if it is enough to test this)
673 if ( !parent->jastdata.supports(parameter_derivatives) ) return 0;
674
675 //cout << "BCS_wf: doing analytic derivatives with respect to parameters."
676 // << endl ;
677
678 int nparms_full=parent->nparms();
679 int nparms_start=0;
680 int nparms_end=nparms_full;
681 int nparms=nparms_end-nparms_start;
682
683 parm_derivatives.gradient.Resize(nparms);
684 parm_derivatives.hessian.Resize(nparms, nparms);
685
686 // gradients of the two-body orbital
687 Array3 <doublevar> twobody_parm_deriv;
688 jast.updateLap(&parent->jastdata,sample);
689 jast.get_twobody_ParmDeriv(sample,twobody_parm_deriv);
690
691 // gradient
692 parm_derivatives.gradient=0.0;
693 for(int det=0; det < ndet; det++) {
694 for (int alpha=0; alpha<nparms; alpha++) {
695 for(int i=0; i< nelectrons(0); i++) {
696 for (int j=0; j< nelectrons(1); j++) {
697 parm_derivatives.gradient(alpha)+=
698 twobody_parm_deriv(i,nelectrons(0)+j,alpha)*inverse(det)(i,j)*
699 parent->magnification_factor;
700 }
701 }
702 }
703 }
704
705 // hessian
706 parm_derivatives.hessian=0.0;
707 for(int det=0; det < ndet; det++) {
708 for (int alpha=0; alpha<nparms; alpha++) {
709 for (int beta=0; beta<nparms; beta++) {
710 parm_derivatives.hessian(alpha,beta)+=
711 parm_derivatives.gradient(alpha)*parm_derivatives.gradient(beta);
712 for(int i=0; i< nelectrons(0); i++) {
713 for (int j=0; j< nelectrons(1); j++) {
714 for(int k=0; k< nelectrons(0); k++) {
715 for (int l=0; l< nelectrons(1); l++) {
716 parm_derivatives.hessian(alpha,beta)-=
717 twobody_parm_deriv(i,nelectrons(0)+j,alpha)*inverse(det)(k,j)*
718 twobody_parm_deriv(k,nelectrons(0)+l,beta)*inverse(det)(i,l)*
719 parent->magnification_factor*parent->magnification_factor;
720 }
721 }
722 }
723 }
724 }
725 }
726 }
727
728 return 1;
729
730 }
731
732
733 //-------------------------------------------------------------------------
734
plot1DInternals(Array1<doublevar> & xdata,vector<Array1<doublevar>> & data,vector<string> & desc,string desc0)735 void BCS_wf::plot1DInternals(Array1 <doublevar> & xdata,
736 vector <Array1 <doublevar> > & data,
737 vector <string> & desc,
738 string desc0) {
739
740 jast.plot1DInternals(xdata,data,desc,"BCS, ");
741
742 }
743
744 //-------------------------------------------------------------------------
745