GPUE  v1.0
GPU Gross-Pitaevskii Equation numerical solver for Bose-Einstein condensates
init.cu
Go to the documentation of this file.
1 #include "../include/init.h"
2 #include "../include/dynamic.h"
3 
4 int init(Grid &par){
5 
6  set_fns(par);
7 
8  // Re-establishing variables from parsed Grid class
9  // Initializes uninitialized variables to 0 values
10  std::string data_dir = par.sval("data_dir");
11  int dimnum = par.ival("dimnum");
12  int N = par.ival("atoms");
13  int xDim = par.ival("xDim");
14  int yDim = par.ival("yDim");
15  int zDim = par.ival("zDim");
16  bool write_file = par.bval("write_file");
17  bool cyl_coord = par.bval("cyl_coord");
18  bool corotating = par.bval("corotating");
19  dim3 threads;
20  unsigned int gSize = xDim;
21  if (dimnum > 1){
22  gSize *= yDim;
23  }
24  if (dimnum > 2){
25  gSize *= zDim;
26  }
27  double omega = par.dval("omega");
28  double gdt = par.dval("gdt");
29  double dt = par.dval("dt");
30  double omegaX = par.dval("omegaX");
31  double omegaY = par.dval("omegaY");
32  double omegaZ = par.dval("omegaZ");
33  double gammaY = par.dval("gammaY"); //Aspect ratio of trapping geometry.
34  double l = par.dval("winding");
35  double box_size = par.dval("box_size");
36  double *Energy;
37  double *r;
38  double *V_opt;
39  double *Energy_gpu;
40  cufftDoubleComplex *wfc;
41  if (par.bval("read_wfc") == true){
42  wfc = par.cufftDoubleComplexval("wfc");
43  }
44  cufftDoubleComplex *EV_opt;
45  cufftDoubleComplex *wfc_backup;
46  cufftDoubleComplex *EappliedField;
47 
48  std::cout << "gSize is: " << gSize << '\n';
49  cufftResult result;
50  cufftHandle plan_1d;
51  cufftHandle plan_2d;
52  cufftHandle plan_3d;
53  cufftHandle plan_other2d;
54  cufftHandle plan_dim2;
55  cufftHandle plan_dim3;
56 
57  std::string buffer;
58  double Rxy; //Condensate scaling factor.
59  double a0x, a0y, a0z; //Harmonic oscillator length in x and y directions
60 
61  generate_grid(par);
62 
63  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
64 
65  double mass = 1.4431607e-25; //Rb 87 mass, kg
66  par.store("mass",mass);
67  double a_s = 4.76e-9;
68  par.store("a_s",a_s);
69 
70  double sum = 0.0;
71 
72  a0x = sqrt(HBAR/(2*mass*omegaX));
73  a0y = sqrt(HBAR/(2*mass*omegaY));
74  a0z = sqrt(HBAR/(2*mass*omegaZ));
75  par.store("a0x",a0x);
76  par.store("a0y",a0y);
77  par.store("a0z",a0z);
78 
79  // Let's go ahead and define the gDensConst here
80  // N*4*HBAR*HBAR*PI*(4.67e-9/mass)*sqrt(mass*(omegaZ)/(2*PI*HBAR)
81  double gDenConst = N*4*HBAR*HBAR*PI*(a_s/mass);
82  if (dimnum == 2){
83  gDenConst*= sqrt(mass*(omegaZ)/(2*PI*HBAR));
84  }
85  par.store("gDenConst", gDenConst);
86 
87  Rxy = pow(15,0.2)*pow(N*a_s*sqrt(mass*omegaZ/HBAR),0.2);
88  par.store("Rxy",Rxy);
89 
90  //std::cout << "Rxy is: " << Rxy << '\n';
91  double xMax, yMax, zMax;
92  if (box_size > 0){
93  xMax = box_size;
94  yMax = box_size;
95  zMax = box_size;
96  }
97  else{
98  xMax = 6*Rxy*a0x;
99  yMax = 6*Rxy*a0y;
100  zMax = 6*Rxy*a0z;
101  }
102  par.store("xMax",xMax);
103  par.store("yMax",yMax);
104  par.store("zMax",zMax);
105 
106  double pxMax, pyMax, pzMax;
107  pxMax = (PI/xMax)*(xDim>>1);
108  pyMax = (PI/yMax)*(yDim>>1);
109  pzMax = (PI/zMax)*(zDim>>1);
110  par.store("pyMax",pyMax);
111  par.store("pxMax",pxMax);
112  par.store("pzMax",pzMax);
113 
114  double dx = xMax/(xDim>>1);
115  double dy = yMax/(yDim>>1);
116  double dz = zMax/(zDim>>1);
117  if (dimnum < 3){
118  dz = 1;
119  }
120  if (dimnum < 2){
121  dy = 1;
122  }
123  par.store("dx",dx);
124  par.store("dy",dy);
125  par.store("dz",dz);
126 
127  double dpx, dpy, dpz;
128  dpx = PI/(xMax);
129  dpy = PI/(yMax);
130  dpz = PI/(zMax);
131  //std::cout << "yMax is: " << yMax << '\t' << "xMax is: " << xMax << '\n';
132  //std::cout << "dpx and dpy are:" << '\n';
133  //std::cout << dpx << '\t' << dpy << '\n';
134  par.store("dpx",dpx);
135  par.store("dpy",dpy);
136  par.store("dpz",dpz);
137 
138 
139  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
140 
141  /* Initialise wavefunction, momentum, position, angular momentum,
142  imaginary and real-time evolution operators . */
143  Energy = (double*) malloc(sizeof(double) * gSize);
144  r = (double *) malloc(sizeof(double) * gSize);
145  V_opt = (double *) malloc(sizeof(double) * gSize);
146  EV_opt = (cufftDoubleComplex *) malloc(sizeof(cufftDoubleComplex) * gSize);
147  EappliedField = (cufftDoubleComplex *) malloc(sizeof(cufftDoubleComplex) *
148  gSize);
149 
150  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
151 
152 /*
153  #ifdef __linux
154  int cores = omp_get_num_procs();
155  par.store("Cores_Total",cores);
156 
157  // Assuming dev system specifics (Xeon with HT -> cores detected / 2)
158  par.store("Cores_Max",cores/2);
159  omp_set_num_threads(cores/2);
160 
161  //#pragma omp parallel for private(j)
162  #endif
163 */
164 
165  par.store("gSize", xDim*yDim*zDim);
166  if (par.bval("use_param_file")){
167  parse_param_file(par);
168  }
169  generate_fields(par);
170  double *K = par.dsval("K");
171  double *Ax = par.dsval("Ax");
172  double *Ay = par.dsval("Ay");
173  double *Az = par.dsval("Az");
174  double *V = par.dsval("V");
175 
176  double *pAx = par.dsval("pAx");
177  double *pAy = par.dsval("pAy");
178  double *pAz = par.dsval("pAz");
179 
180  double *x = par.dsval("x");
181  double *y = par.dsval("y");
182  double *z = par.dsval("z");
183 
184  double2 *GpAx = par.cufftDoubleComplexval("GpAx");
185  double2 *GpAy = par.cufftDoubleComplexval("GpAy");
186  double2 *GpAz = par.cufftDoubleComplexval("GpAz");
187  double2 *EpAx = par.cufftDoubleComplexval("EpAx");
188  double2 *EpAy = par.cufftDoubleComplexval("EpAy");
189  double2 *EpAz = par.cufftDoubleComplexval("EpAz");
190 
191  double2 *GV = par.cufftDoubleComplexval("GV");
192  double2 *EV = par.cufftDoubleComplexval("EV");
193  double2 *GK = par.cufftDoubleComplexval("GK");
194  double2 *EK = par.cufftDoubleComplexval("EK");
195 
196  wfc = par.cufftDoubleComplexval("wfc");
197 
198  int index = 0;
199  for(int i=0; i < gSize; i++ ){
200  sum+=sqrt(wfc[i].x*wfc[i].x + wfc[i].y*wfc[i].y);
201  }
202 
203  if (write_file){
204  double *Bz;
205  double *Bx;
206  double *By;
207  if (dimnum == 2){
208  Bz = curl2d(par, Ax, Ay);
209  }
210  if (dimnum == 3){
211  std::cout << "Calculating the 3d curl..." << '\n';
212  Bx = curl3d_x(par, Ax, Ay, Az);
213  By = curl3d_y(par, Ax, Ay, Az);
214  Bz = curl3d_z(par, Ax, Ay, Az);
215  std::cout << "Finished calculating Curl" << '\n';
216  }
217  std::cout << "writing initial variables to file..." << '\n';
218  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
219  //hdfWriteDouble(xDim, V, 0, "V_0"); //HDF COMING SOON!
220  //hdfWriteComplex(xDim, wfc, 0, "wfc_0");
221  if (cyl_coord && dimnum > 2){
222  double *Br = curl3d_r(par, Bx, By);
223  double *Bphi = curl3d_phi(par, Bx, By);
224 
225  FileIO::writeOutDouble(buffer, data_dir + "Br",Br,gSize,0);
226  FileIO::writeOutDouble(buffer, data_dir + "Bphi",Bphi,gSize,0);
227  FileIO::writeOutDouble(buffer, data_dir + "Bz",Bz,gSize,0);
228 
229  free(Br);
230  free(Bx);
231  free(By);
232  free(Bz);
233  free(Bphi);
234  }
235  else{
236  if (dimnum > 1){
237  FileIO::writeOutDouble(buffer, data_dir + "Bz",Bz,gSize,0);
238  free(Bz);
239  }
240  if (dimnum > 2){
241  FileIO::writeOutDouble(buffer, data_dir + "Bx",Bx,gSize,0);
242  FileIO::writeOutDouble(buffer, data_dir + "By",By,gSize,0);
243  free(Bx);
244  free(By);
245  }
246  }
247 
248  FileIO::writeOutDouble(buffer, data_dir + "V",V,gSize,0);
249  FileIO::writeOutDouble(buffer, data_dir + "K",K,gSize,0);
250  FileIO::writeOutDouble(buffer, data_dir + "pAy",pAy,gSize,0);
251  FileIO::writeOutDouble(buffer, data_dir + "pAx",pAx,gSize,0);
252  FileIO::writeOutDouble(buffer, data_dir + "Ax",Ax,gSize,0);
253  FileIO::writeOutDouble(buffer, data_dir + "Ay",Ay,gSize,0);
254  FileIO::writeOutDouble(buffer, data_dir + "Az",Az,gSize,0);
255  FileIO::writeOutDouble(buffer, data_dir + "x",x,xDim,0);
256  FileIO::writeOutDouble(buffer, data_dir + "y",y,yDim,0);
257  FileIO::writeOutDouble(buffer, data_dir + "z",z,zDim,0);
258  FileIO::writeOut(buffer, data_dir + "WFC",wfc,gSize,0);
259  FileIO::writeOut(buffer, data_dir + "EpAz",EpAz,gSize,0);
260  FileIO::writeOut(buffer, data_dir + "EpAy",EpAy,gSize,0);
261  FileIO::writeOut(buffer, data_dir + "EpAx",EpAx,gSize,0);
262  FileIO::writeOut(buffer, data_dir + "GK",GK,gSize,0);
263  FileIO::writeOut(buffer, data_dir + "GV",GV,gSize,0);
264  FileIO::writeOut(buffer, data_dir + "GpAx",GpAx,gSize,0);
265  FileIO::writeOut(buffer, data_dir + "GpAy",GpAy,gSize,0);
266  FileIO::writeOut(buffer, data_dir + "GpAz",GpAz,gSize,0);
267  }
268 
269  if (par.bval("read_wfc") == false){
270  sum=sqrt(sum*dx*dy*dz);
271  for (int i = 0; i < gSize; i++){
272  wfc[i].x = (wfc[i].x)/(sum);
273  wfc[i].y = (wfc[i].y)/(sum);
274  }
275  }
276 
277  result = cufftPlan2d(&plan_2d, xDim, yDim, CUFFT_Z2Z);
278  if(result != CUFFT_SUCCESS){
279  printf("Result:=%d\n",result);
280  printf("Error: Could not execute cufftPlan2d(%s, %d, %d).\n", "plan_2d",
281  (unsigned int)xDim, (unsigned int)yDim);
282  exit(1);
283  }
284  generate_plan_other2d(&plan_other2d, par);
285 
286  generate_plan_other3d(&plan_1d, par, 0);
287  generate_plan_other3d(&plan_dim2, par, 1);
288  generate_plan_other3d(&plan_dim3, par, 2);
289  result = cufftPlan3d(&plan_3d, xDim, yDim, zDim, CUFFT_Z2Z);
290  if(result != CUFFT_SUCCESS){
291  printf("Result:=%d\n",result);
292  printf("Error: Could not execute cufftPlan3d(%s, %d, %d, %d).\n",
293  "plan_3d",
294  (unsigned int)xDim, (unsigned int)yDim, (unsigned int) zDim);
295  exit(1);
296  }
297 
298  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
299 
300  //std::cout << GV[0].x << '\t' << GK[0].x << '\t'
301  // << pAy[0] << '\t' << pAx[0] << '\n';
302 
303  //std::cout << "storing variables..." << '\n';
304 
305  // Storing variables that have been initialized
306  // Re-establishing variables from parsed Grid class
307  // Initializes uninitialized variables to 0 values
308  par.store("Energy", Energy);
309  par.store("r", r);
310  par.store("Energy_gpu", Energy_gpu);
311  par.store("wfc", wfc);
312  par.store("EV_opt", EV_opt);
313  par.store("V_opt", V_opt);
314  par.store("wfc_backup", wfc_backup);
315  par.store("EappliedField", EappliedField);
316 
317  par.store("result", result);
318  par.store("plan_1d", plan_1d);
319  par.store("plan_2d", plan_2d);
320  par.store("plan_other2d", plan_other2d);
321  par.store("plan_3d", plan_3d);
322  par.store("plan_dim2", plan_dim2);
323  par.store("plan_dim3", plan_dim3);
324 
325  // Parameters for time-depd variables.
326  par.store("K_time", false);
327  par.store("V_time", false);
328  par.store("Ax_time", false);
329  par.store("Ay_time", false);
330  par.store("Az_time", false);
331 
332 
333  std::cout << "variables stored" << '\n';
334 
335  return 0;
336 }
337 
338 void set_variables(Grid &par, bool ev_type){
339  // Re-establishing variables from parsed Grid class
340  // Note that 3d variables are set to nullptr's unless needed
341  // This might need to be fixed later
342  double dx = par.dval("dx");
343  double dy = par.dval("dy");
344  double *V_opt = par.dsval("V_opt");
345  double *pAy = par.dsval("pAy");
346  double *pAx = par.dsval("pAx");
347  double2 *pAy_gpu;
348  double2 *pAx_gpu;
349  double2 *pAz_gpu;
350  double2 *V_gpu;
351  double2 *K_gpu;
352  cufftDoubleComplex *wfc = par.cufftDoubleComplexval("wfc");
353  cufftDoubleComplex *wfc_gpu = par.cufftDoubleComplexval("wfc_gpu");
354  cudaError_t err;
355  int dimnum = par.ival("dimnum");
356  int xDim = par.ival("xDim");
357  int yDim = par.ival("yDim");
358  int zDim = par.ival("zDim");
359  int gsize = xDim;
360 
361  // Special variables for the 3d case
362  if (dimnum > 1){
363  gsize *= yDim;
364  }
365  if (dimnum > 2){
366  gsize *= zDim;
367  }
368  if(!par.bval("V_time")){
369  cudaMalloc((void**) &V_gpu, sizeof(double2)*gsize);
370  }
371  if(!par.bval("K_time")){
372  cudaMalloc((void**) &K_gpu, sizeof(double2)*gsize);
373  }
374  if(!par.bval("Ax_time")){
375  cudaMalloc((void**) &pAx_gpu, sizeof(double2)*gsize);
376  }
377  if(!par.bval("Ay_time") && dimnum > 1){
378  cudaMalloc((void**) &pAy_gpu, sizeof(double2)*gsize);
379  }
380  if(!par.bval("Az_time") && dimnum > 2){
381  cudaMalloc((void**) &pAz_gpu, sizeof(double2)*gsize);
382  }
383 
384  if (ev_type == 0){
385  cufftDoubleComplex *GK = par.cufftDoubleComplexval("GK");
386  cufftDoubleComplex *GV = par.cufftDoubleComplexval("GV");
387  cufftDoubleComplex *GpAx = par.cufftDoubleComplexval("GpAx");
388  cufftDoubleComplex *GpAy = nullptr;
389  cufftDoubleComplex *GpAz = nullptr;
390 
391  if(!par.bval("K_time")){
392  err=cudaMemcpy(K_gpu, GK, sizeof(cufftDoubleComplex)*gsize,
393  cudaMemcpyHostToDevice);
394  if(err!=cudaSuccess){
395  std::cout << "ERROR: Could not copy K_gpu to device" << '\n';
396  exit(1);
397  }
398  }
399  if(!par.bval("V_time")){
400  err=cudaMemcpy(V_gpu, GV, sizeof(cufftDoubleComplex)*gsize,
401  cudaMemcpyHostToDevice);
402  if(err!=cudaSuccess){
403  std::cout << "ERROR: Could not copy V_gpu to device" << '\n';
404  exit(1);
405  }
406  }
407  if(!par.bval("Ax_time")){
408  err=cudaMemcpy(pAx_gpu, GpAx, sizeof(cufftDoubleComplex)*gsize,
409  cudaMemcpyHostToDevice);
410  if(err!=cudaSuccess){
411  std::cout << "ERROR: Could not copy pAx_gpu to device" << '\n';
412  exit(1);
413  }
414  }
415  err=cudaMemcpy(wfc_gpu, wfc, sizeof(cufftDoubleComplex)*gsize,
416  cudaMemcpyHostToDevice);
417  if(err!=cudaSuccess){
418  std::cout << "ERROR: Could not copy wfc_gpu to device" << '\n';
419  exit(1);
420  }
421  par.store("K_gpu", K_gpu);
422  par.store("V_gpu", V_gpu);
423  par.store("wfc_gpu", wfc_gpu);
424  par.store("pAy_gpu", pAy_gpu);
425  par.store("pAx_gpu", pAx_gpu);
426 
427  // Special cases for 3d
428  if (dimnum > 1 && !par.bval("Ay_time")){
429  GpAy = par.cufftDoubleComplexval("GpAy");
430  err=cudaMemcpy(pAy_gpu, GpAy, sizeof(cufftDoubleComplex)*gsize,
431  cudaMemcpyHostToDevice);
432 
433  if(err!=cudaSuccess){
434  std::cout << "ERROR: Could not copy pAy_gpu to device" << '\n';
435  exit(1);
436  }
437  par.store("pAy_gpu", pAy_gpu);
438 
439  }
440  if (dimnum > 2 && !par.bval("Az_time")){
441  GpAz = par.cufftDoubleComplexval("GpAz");
442  err=cudaMemcpy(pAz_gpu, GpAz, sizeof(cufftDoubleComplex)*gsize,
443  cudaMemcpyHostToDevice);
444 
445  if(err!=cudaSuccess){
446  std::cout << "ERROR: Could not copy pAz_gpu to device" << '\n';
447  exit(1);
448  }
449  par.store("pAz_gpu", pAz_gpu);
450 
451  }
452  free(GV); free(GK); free(GpAy); free(GpAx); free(GpAz);
453  }
454  else if (ev_type == 1){
455 
456  cufftDoubleComplex *EV = par.cufftDoubleComplexval("EV");
457  cufftDoubleComplex *EK = par.cufftDoubleComplexval("EK");
458  cufftDoubleComplex *EpAx = par.cufftDoubleComplexval("EpAx");
459  cufftDoubleComplex *EpAy = nullptr;
460  cufftDoubleComplex *EpAz = nullptr;
461  if (!par.bval("K_time")){
462  err=cudaMemcpy(K_gpu, EK, sizeof(cufftDoubleComplex)*gsize,
463  cudaMemcpyHostToDevice);
464  if(err!=cudaSuccess){
465  std::cout << "ERROR: Could not copy K_gpu to device" << '\n';
466  exit(1);
467  }
468  par.store("K_gpu", K_gpu);
469  }
470  if(!par.bval("Ax_time")){
471  err=cudaMemcpy(pAx_gpu, EpAx, sizeof(cufftDoubleComplex)*gsize,
472  cudaMemcpyHostToDevice);
473  if(err!=cudaSuccess){
474  std::cout << "ERROR: Could not copy pAx_gpu to device" << '\n';
475  std::cout << err << '\n';
476  exit(1);
477  }
478  par.store("pAx_gpu", pAx_gpu);
479  }
480 
481  if (!par.bval("V_time")){
482  err=cudaMemcpy(V_gpu, EV, sizeof(cufftDoubleComplex)*gsize,
483  cudaMemcpyHostToDevice);
484  if(err!=cudaSuccess){
485  std::cout << "ERROR: Could not copy V_gpu to device" << '\n';
486  exit(1);
487  }
488  par.store("V_gpu", V_gpu);
489  }
490  err=cudaMemcpy(wfc_gpu, wfc, sizeof(cufftDoubleComplex)*gsize,
491  cudaMemcpyHostToDevice);
492  if(err!=cudaSuccess){
493  std::cout << "ERROR: Could not copy wfc_gpu to device" << '\n';
494  exit(1);
495  }
496 
497  par.store("wfc_gpu", wfc_gpu);
498 
499  // Special variables / instructions for 2/3d case
500  if (dimnum > 1 && !par.bval("Ay_time")){
501  pAy_gpu = par.cufftDoubleComplexval("pAy_gpu");
502  EpAy = par.cufftDoubleComplexval("EpAy");
503  err=cudaMemcpy(pAy_gpu, EpAy, sizeof(cufftDoubleComplex)*gsize,
504  cudaMemcpyHostToDevice);
505  if(err!=cudaSuccess){
506  std::cout << "ERROR: Could not copy pAy_gpu to device" << '\n';
507  exit(1);
508  }
509  par.store("pAy_gpu", pAy_gpu);
510  }
511 
512  if (dimnum > 2 && !par.bval("Az_time")){
513  pAz_gpu = par.cufftDoubleComplexval("pAz_gpu");
514  EpAz = par.cufftDoubleComplexval("EpAz");
515  err=cudaMemcpy(pAz_gpu, EpAz, sizeof(cufftDoubleComplex)*gsize,
516  cudaMemcpyHostToDevice);
517  if(err!=cudaSuccess){
518  std::cout << "ERROR: Could not copy pAz_gpu to device" << '\n';
519  exit(1);
520  }
521  par.store("pAz_gpu", pAz_gpu);
522  }
523 
524  free(EV); free(EK); free(EpAy); free(EpAx); free(EpAz);
525  }
526 
527 }
528 
529 int main(int argc, char **argv){
530 
531  Grid par = parseArgs(argc,argv);
532  //Grid par2 = parseArgs(argc,argv);
533 
534  int device = par.ival("device");
535  int dimnum = par.ival("dimnum");
536  cudaSetDevice(device);
537 
538  std::string buffer;
539  time_t start,fin;
540  time(&start);
541  printf("Start: %s\n", ctime(&start));
542 
543  //************************************************************//
544  /*
545  * Initialise the Params data structure to track params and variables
546  */
547  //************************************************************//
548 
549  // If we want to read in a wfc, we may also need to imprint a phase. This
550  // will be done in the init_2d and init_3d functions
551  // We need a number of parameters for now
552  int xDim = par.ival("xDim");
553  int yDim = par.ival("yDim");
554  int zDim = par.ival("zDim");
555  if(par.bval("read_wfc") == true){
556 
557  // Initializing the wfc
558  int gSize = xDim * yDim * zDim;
559  cufftDoubleComplex *wfc;
560 
561  std::string infile = par.sval("infile");
562  std::string infilei = par.sval("infilei");
563  printf("Loading wavefunction...");
564  wfc=FileIO::readIn(infile,infilei,gSize);
565  par.store("wfc",wfc);
566  printf("Wavefunction loaded.\n");
567  //std::string data_dir = par.sval("data_dir");
568  //FileIO::writeOut(buffer, data_dir + "WFC_CHECK",wfc,gSize,0);
569  }
570 
571  init(par);
572 
573  int gsteps = par.ival("gsteps");
574  int esteps = par.ival("esteps");
575  std::string data_dir = par.sval("data_dir");
576  std::cout << "variables re-established" << '\n';
577 
578  if (par.bval("write_file")){
579  FileIO::writeOutParam(buffer, par, data_dir + "Params.dat");
580  }
581 
582  if(gsteps > 0){
583  std::cout << "Imaginary-time evolution started..." << '\n';
584  set_variables(par, 0);
585 
586  evolve(par, gsteps, 0, buffer);
587  }
588 
589  if(esteps > 0){
590  std::cout << "real-time evolution started..." << '\n';
591  set_variables(par, 1);
592  evolve(par, esteps, 1, buffer);
593  }
594 
595  std::cout << "done evolving" << '\n';
596  time(&fin);
597  printf("Finish: %s\n", ctime(&fin));
598  printf("Total time: %ld seconds\n ",(long)fin-start);
599  std::cout << '\n';
600  return 0;
601 }
602