Changeset 6edeca for pcp/src


Ignore:
Timestamp:
Apr 22, 2008, 9:14:25 AM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
b0aa9c
Parents:
0f2c43
Message:

Moved "testing" DoBrent minimisation of perturbed psis to end of file

Functions PerturbedFunction_..f..() were shifted to end of file, MinimisePerturbed() was copied and put into comment block. In the original DoBrent lines (all GSL minimisation stuff) was removed.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified pcp/src/perturbed.c

    r0f2c43 r6edeca  
    6969#include "wannier.h"
    7070
    71 /** evaluates perturbed energy functional.
    72  * \param norm norm of current Psi in functional
    73  * \param *params void-pointer to parameter array
    74  * \return evaluated functional at f(x) with \a norm
    75  */
    76 double perturbed_function (double norm, void *params) {
    77   struct Problem *P = (struct Problem *)params;
    78   int i, n = P->R.LevS->MaxG;
    79   double old_norm  = GramSchGetNorm2(P,P->R.LevS,P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo]);
    80   fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
    81   fprintf(stderr,"(%i) perturbed_function: setting norm to %lg ...", P->Par.me, norm);
    82   // set desired norm for current Psi
    83   for (i=0; i< n; i++) {
    84     currentPsi[i].re *= norm/old_norm; // real part
    85     currentPsi[i].im *= norm/old_norm; // imaginary part
    86   }
    87   P->R.PsiStep = 0; // make it not advance to next Psi
    88  
    89   //debug(P,"UpdateActualPsiNo");
    90   UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
    91   //debug(P,"UpdateEnergyArray");
    92   UpdateEnergyArray(P); // shift energy values in their array by one
    93   //debug(P,"UpdatePerturbedEnergyCalculation");
    94   UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
    95   EnergyAllReduce(P); // gather from all processes and sum up to total energy
    96 /*
    97   for (i=0; i< n; i++) {
    98     currentPsi[i].re /= norm/old_norm; // real part
    99     currentPsi[i].im /= norm/old_norm; // imaginary part
    100   }*/
    101  
    102   fprintf(stderr,"%lg\n", P->Lat.E->TotalEnergy[0]);
    103   return P->Lat.E->TotalEnergy[0];   // and return evaluated functional
    104 }
    105 
    106 /** evaluates perturbed energy functional.
    107  * \param *x current position in functional
    108  * \param *params void-pointer to parameter array
    109  * \return evaluated functional at f(x)
    110  */
    111 double perturbed_f (const gsl_vector *x, void *params) {
    112   struct Problem *P = (struct Problem *)params;
    113   int i, n = P->R.LevS->MaxG*2;
    114   fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
    115   //int diff = 0;
    116   //debug(P,"f");
    117   // put x into current Psi
    118   for (i=0; i< n; i+=2) {
    119     //if ((currentPsi[i/2].re != gsl_vector_get (x, i)) || (currentPsi[i/2].im != gsl_vector_get (x, i+1))) diff++;
    120     currentPsi[i/2].re = gsl_vector_get (x, i);   // real part
    121     currentPsi[i/2].im = gsl_vector_get (x, i+1); // imaginary part
    122   }
    123   //if (diff) fprintf(stderr,"(%i) %i differences between old and new currentPsi.\n", P->Par.me, diff);
    124   P->R.PsiStep = 0; // make it not advance to next Psi
    125  
    126   //debug(P,"UpdateActualPsiNo");
    127   UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
    128   //debug(P,"UpdateEnergyArray");
    129   UpdateEnergyArray(P); // shift energy values in their array by one
    130   //debug(P,"UpdatePerturbedEnergyCalculation");
    131   UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
    132   EnergyAllReduce(P); // gather from all processes and sum up to total energy
    133  
    134   return P->Lat.E->TotalEnergy[0];   // and return evaluated functional
    135 }
    136 
    137 /** evaluates perturbed energy gradient.
    138  * \param *x current position in functional
    139  * \param *params void-pointer to parameter array
    140  * \param *g array for gradient vector on return
    141  */
    142 void perturbed_df (const gsl_vector *x, void *params, gsl_vector *g) {
    143   struct Problem *P = (struct Problem *)params;
    144   int i, n = P->R.LevS->MaxG*2;
    145   fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
    146   fftw_complex *gradient = P->Grad.GradientArray[ActualGradient];
    147   //int diff = 0;
    148   //debug(P,"df");
    149   // put x into current Psi
    150   for (i=0; i< n; i+=2) {
    151     //if ((currentPsi[i/2].re != gsl_vector_get (x, i)) || (currentPsi[i/2].im != gsl_vector_get (x, i+1))) diff++;
    152     currentPsi[i/2].re = gsl_vector_get (x, i);   // real part
    153     currentPsi[i/2].im = gsl_vector_get (x, i+1); // imaginary part
    154   }
    155   //if (diff) fprintf(stderr,"(%i) %i differences between old and new currentPsi.\n", P->Par.me, diff);
    156   P->R.PsiStep = 0; // make it not advance to next Psi
    157  
    158   //debug(P,"UpdateActualPsiNo");
    159   UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
    160   //debug(P,"UpdateEnergyArray");
    161   UpdateEnergyArray(P); // shift energy values in their array by one
    162   //debug(P,"UpdatePerturbedEnergyCalculation");
    163   UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
    164   EnergyAllReduce(P); // gather from all processes and sum up to total energy
    165  
    166   // checkout gradient
    167   //diff = 0;
    168   for (i=0; i< n; i+=2) {
    169     //if ((-gradient[i/2].re != gsl_vector_get (g, i)) || (-gradient[i/2].im != gsl_vector_get (g, i+1)))  diff++;
    170     gsl_vector_set (g, i, -gradient[i/2].re);   // real part
    171     gsl_vector_set (g, i+1, -gradient[i/2].im); // imaginary part
    172   }
    173   //if (diff) fprintf(stderr,"(%i) %i differences between old and new gradient.\n", P->Par.me, diff);
    174 }
    175 
    176 /** evaluates perturbed energy functional and gradient.
    177  * \param *x current position in functional
    178  * \param *params void-pointer to parameter array
    179  * \param *f pointer to energy function value on return
    180  * \param *g array for gradient vector on return
    181  */
    182 void perturbed_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *g) {
    183   struct Problem *P = (struct Problem *)params;
    184   int i, n = P->R.LevS->MaxG*2;
    185   fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
    186   fftw_complex *gradient = P->Grad.GradientArray[ActualGradient];
    187   //int diff = 0;
    188   //debug(P,"fdf");
    189   // put x into current Psi
    190   for (i=0; i< n; i+=2) {
    191     //if ((currentPsi[i/2].re != gsl_vector_get (x, i)) || (currentPsi[i/2].im != gsl_vector_get (x, i+1))) diff++;
    192     currentPsi[i/2].re = gsl_vector_get (x, i);   // real part
    193     currentPsi[i/2].im = gsl_vector_get (x, i+1); // imaginary part
    194   }
    195   //if (diff) fprintf(stderr,"(%i) %i differences between old and new currentPsi.\n", P->Par.me, diff);
    196   P->R.PsiStep = 0; // make it not advance to next Psi
    197  
    198   //debug(P,"UpdateActualPsiNo");
    199   UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
    200   //debug(P,"UpdateEnergyArray");
    201   UpdateEnergyArray(P); // shift energy values in their array by one
    202   //debug(P,"UpdatePerturbedEnergyCalculation");
    203   UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
    204   EnergyAllReduce(P); // gather from all processes and sum up to total energy
    205  
    206   // checkout gradient
    207   //diff = 0;
    208   for (i=0; i< n; i+=2) {
    209     //if ((-gradient[i/2].re != gsl_vector_get (g, i)) || (-gradient[i/2].im != gsl_vector_get (g, i+1)))  diff++;
    210     gsl_vector_set (g, i, -gradient[i/2].re);   // real part
    211     gsl_vector_set (g, i+1, -gradient[i/2].im); // imaginary part
    212   }
    213   //if (diff) fprintf(stderr,"(%i) %i differences between old and new gradient.\n", P->Par.me, diff);
    214  
    215   *f = P->Lat.E->TotalEnergy[0];   // and return evaluated functional
    216 }
    21771
    21872/** Minimisation of the PsiTypeTag#Perturbed_RxP0, PsiTypeTag#Perturbed_P0 and other orbitals.
    … …  
    258112  struct Psis *Psi = &Lat->Psi;
    259113  int type;
    260   //int i;
    261  
    262   // stuff for GSL minimization
    263   //size_t iter;
    264   //int status, Status
    265   int n = R->LevS->MaxG*2;
    266   const gsl_multimin_fdfminimizer_type *T_multi;
    267   const gsl_min_fminimizer_type *T;
    268   gsl_multimin_fdfminimizer *s_multi;
    269   gsl_min_fminimizer *s;
    270   gsl_vector *x;//, *ss;
    271   gsl_multimin_function_fdf my_func;
    272   gsl_function F;
    273   //fftw_complex *currentPsi;
    274   //double a,b,m, f_m, f_a, f_b;
    275   //double old_norm;
    276  
    277   my_func.f = &perturbed_f;
    278   my_func.df = &perturbed_df;
    279   my_func.fdf = &perturbed_fdf;
    280   my_func.n = n;
    281   my_func.params = P;
    282   F.function = &perturbed_function;
    283   F.params = P;
    284      
    285   x = gsl_vector_alloc (n);
    286   //ss = gsl_vector_alloc (Psi->NoOfPsis);
    287   T_multi = gsl_multimin_fdfminimizer_vector_bfgs;
    288   s_multi = gsl_multimin_fdfminimizer_alloc (T_multi, n);
    289   T = gsl_min_fminimizer_brent;
    290   s = gsl_min_fminimizer_alloc (T);
    291114 
    292115  for (type=Perturbed_P0;type<=Perturbed_RxP2;type++) {  // go through each perturbation group separately //
    … …  
    339162      EnergyOutput(P,0);
    340163      while (*Stop != 1) {
    341 /*        // copy current Psi into starting vector
    342         currentPsi = R->LevS->LPsi->LocalPsi[R->ActualLocalPsiNo];
    343         for (i=0; i< n; i+=2) {
    344           gsl_vector_set (x, i, currentPsi[i/2].re);   // real part
    345           gsl_vector_set (x, i+1, currentPsi[i/2].im); // imaginary part
    346         }
    347         gsl_multimin_fdfminimizer_set (s_multi, &my_func, x, 0.01, 1e-2); 
    348         iter = 0;
    349         status = 0;
    350         do {  // look for minimum along current local psi
    351           iter++;
    352           status = gsl_multimin_fdfminimizer_iterate (s_multi);
    353           MPI_Allreduce(&status, &Status, 1, MPI_INT, MPI_MAX, P->Par.comm_ST_Psi);
    354           if (Status)
    355             break;
    356           status = gsl_multimin_test_gradient (s_multi->gradient, 1e-2);
    357           MPI_Allreduce(&status, &Status, 1, MPI_INT, MPI_MAX, P->Par.comm_ST_Psi);
    358           //if (Status == GSL_SUCCESS)
    359             //printf ("Minimum found at:\n");
    360           if (P->Par.me == 0) fprintf (stderr,"(%i,%i,%i)S(%i,%i,%i):\t %5d %10.5f\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, (int)iter, s_multi->f);
    361           //TestGramSch(P,R->LevS,Psi, type); // functions are orthonormal?
    362         } while (Status == GSL_CONTINUE && iter < 3);
    363 /*        // now minimize norm of currentPsi (one-dim)
    364         if (0) {
    365           iter = 0;
    366         status = 0;
    367         m = 1.;
    368         a = MYEPSILON;
    369         b = 100.;
    370         f_a = perturbed_function (a, P);
    371         f_b = perturbed_function (b, P);
    372         f_m = perturbed_function (m, P);
    373         //if ((f_m < f_a) && (f_m < f_b)) {
    374           gsl_min_fminimizer_set (s, &F, m, a, b);
    375           do {  // look for minimum along current local psi
    376             iter++;
    377             status = gsl_min_fminimizer_iterate (s);
    378             m = gsl_min_fminimizer_x_minimum (s);
    379             a = gsl_min_fminimizer_x_lower (s);
    380             b = gsl_min_fminimizer_x_upper (s);         
    381             status = gsl_min_test_interval (a, b, 0.001, 0.0);
    382             if (status == GSL_SUCCESS)
    383               printf ("Minimum found at:\n");
    384             printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
    385                (int) iter, a, b,
    386                m, b - a);
    387           } while (status == GSL_CONTINUE && iter < 100);
    388           old_norm = GramSchGetNorm2(P,P->R.LevS,P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo]);
    389           for (i=0; i< n; i++) {
    390             currentPsi[i].re *= m/old_norm; // real part
    391             currentPsi[i].im *= m/old_norm; // imaginary part
    392           }
    393         } else debug(P,"Norm not minimizable!");*/
    394164        //P->R.PsiStep = P->R.MaxPsiStep; // make it advance to next Psi
    395165        FindPerturbedMinimum(P);       
    … …  
    410180      }
    411181      // now release normalization condition and minimize wrt to norm
    412  /*     *Stop = 0;
    413       while (*Stop != 1) {
    414         currentPsi = R->LevS->LPsi->LocalPsi[R->ActualLocalPsiNo];
    415         iter = 0;
    416         status = 0;
    417         m = 1.;
    418         a = 0.001;
    419         b = 10.;
    420         f_a = perturbed_function (a, P);
    421         f_b = perturbed_function (b, P);
    422         f_m = perturbed_function (m, P);
    423         if ((f_m < f_a) && (f_m < f_b)) {
    424           gsl_min_fminimizer_set (s, &F, m, a, b);
    425           do {  // look for minimum along current local psi
    426             iter++;
    427             status = gsl_min_fminimizer_iterate (s);
    428             m = gsl_min_fminimizer_x_minimum (s);
    429             a = gsl_min_fminimizer_x_lower (s);
    430             b = gsl_min_fminimizer_x_upper (s);         
    431             status = gsl_min_test_interval (a, b, 0.001, 0.0);
    432             if (status == GSL_SUCCESS)
    433               printf ("Minimum found at:\n");
    434             printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
    435                (int) iter, a, b,
    436                m, b - a);
    437           } while (status == GSL_CONTINUE && iter < 100);
    438           old_norm = GramSchGetNorm2(P,P->R.LevS,P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo]);
    439           for (i=0; i< n; i++) {
    440             currentPsi[i].re *= m/old_norm; // real part
    441             currentPsi[i].im *= m/old_norm; // imaginary part
    442           }
    443         }
    444         P->R.PsiStep = P->R.MaxPsiStep; // make it advance to next Psi
    445         //debug(P,"UpdateActualPsiNo");
    446         UpdateActualPsiNo(P, type); // step on to next perturbed Psi
    447         if (*SuperStop != 1)
    448           *SuperStop = CheckCPULIM(P);
    449         *Stop = CalculateMinimumStop(P, *SuperStop);
    450         P->Speed.Steps++; // step on
    451         R->LevS->Step++;
    452       }*/
    453182      if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Write %s srcpsi to disk\n", P->Par.me, R->MinimisationName[type]);
    454183      OutputSrcPsiDensity(P, type);
    … …  
    477206  }
    478207  UpdateActualPsiNo(P, Occupied); // step on back to an occupied one
    479  
    480   gsl_multimin_fdfminimizer_free (s_multi);
    481   gsl_min_fminimizer_free (s);
    482   gsl_vector_free (x);
    483   //gsl_vector_free (ss);
    484208}
    485209
    … …  
    39363660  // finished.
    39373661}
     3662
     3663/** evaluates perturbed energy functional
     3664 * \param norm norm of current Psi in functional
     3665 * \param *params void-pointer to parameter array
     3666 * \return evaluated functional at f(x) with \a norm
     3667 */
     3668double perturbed_function (double norm, void *params) {
     3669  struct Problem *P = (struct Problem *)params;
     3670  int i, n = P->R.LevS->MaxG;
     3671  double old_norm  = GramSchGetNorm2(P,P->R.LevS,P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo]);
     3672  fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
     3673  fprintf(stderr,"(%i) perturbed_function: setting norm to %lg ...", P->Par.me, norm);
     3674  // set desired norm for current Psi
     3675  for (i=0; i< n; i++) {
     3676    currentPsi[i].re *= norm/old_norm; // real part
     3677    currentPsi[i].im *= norm/old_norm; // imaginary part
     3678  }
     3679  P->R.PsiStep = 0; // make it not advance to next Psi
     3680 
     3681  //debug(P,"UpdateActualPsiNo");
     3682  UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
     3683  //debug(P,"UpdateEnergyArray");
     3684  UpdateEnergyArray(P); // shift energy values in their array by one
     3685  //debug(P,"UpdatePerturbedEnergyCalculation");
     3686  UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
     3687  EnergyAllReduce(P); // gather from all processes and sum up to total energy
     3688/*
     3689  for (i=0; i< n; i++) {
     3690    currentPsi[i].re /= norm/old_norm; // real part
     3691    currentPsi[i].im /= norm/old_norm; // imaginary part
     3692  }*/
     3693 
     3694  fprintf(stderr,"%lg\n", P->Lat.E->TotalEnergy[0]);
     3695  return P->Lat.E->TotalEnergy[0];   // and return evaluated functional
     3696}
     3697
     3698/** evaluates perturbed energy functional.
     3699 * \param *x current position in functional
     3700 * \param *params void-pointer to parameter array
     3701 * \return evaluated functional at f(x)
     3702 */
     3703double perturbed_f (const gsl_vector *x, void *params) {
     3704  struct Problem *P = (struct Problem *)params;
     3705  int i, n = P->R.LevS->MaxG*2;
     3706  fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
     3707  //int diff = 0;
     3708  //debug(P,"f");
     3709  // put x into current Psi
     3710  for (i=0; i< n; i+=2) {
     3711    //if ((currentPsi[i/2].re != gsl_vector_get (x, i)) || (currentPsi[i/2].im != gsl_vector_get (x, i+1))) diff++;
     3712    currentPsi[i/2].re = gsl_vector_get (x, i);   // real part
     3713    currentPsi[i/2].im = gsl_vector_get (x, i+1); // imaginary part
     3714  }
     3715  //if (diff) fprintf(stderr,"(%i) %i differences between old and new currentPsi.\n", P->Par.me, diff);
     3716  P->R.PsiStep = 0; // make it not advance to next Psi
     3717 
     3718  //debug(P,"UpdateActualPsiNo");
     3719  UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
     3720  //debug(P,"UpdateEnergyArray");
     3721  UpdateEnergyArray(P); // shift energy values in their array by one
     3722  //debug(P,"UpdatePerturbedEnergyCalculation");
     3723  UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
     3724  EnergyAllReduce(P); // gather from all processes and sum up to total energy
     3725 
     3726  return P->Lat.E->TotalEnergy[0];   // and return evaluated functional
     3727}
     3728
     3729/** evaluates perturbed energy gradient.
     3730 * \param *x current position in functional
     3731 * \param *params void-pointer to parameter array
     3732 * \param *g array for gradient vector on return
     3733 */
     3734void perturbed_df (const gsl_vector *x, void *params, gsl_vector *g) {
     3735  struct Problem *P = (struct Problem *)params;
     3736  int i, n = P->R.LevS->MaxG*2;
     3737  fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
     3738  fftw_complex *gradient = P->Grad.GradientArray[ActualGradient];
     3739  //int diff = 0;
     3740  //debug(P,"df");
     3741  // put x into current Psi
     3742  for (i=0; i< n; i+=2) {
     3743    //if ((currentPsi[i/2].re != gsl_vector_get (x, i)) || (currentPsi[i/2].im != gsl_vector_get (x, i+1))) diff++;
     3744    currentPsi[i/2].re = gsl_vector_get (x, i);   // real part
     3745    currentPsi[i/2].im = gsl_vector_get (x, i+1); // imaginary part
     3746  }
     3747  //if (diff) fprintf(stderr,"(%i) %i differences between old and new currentPsi.\n", P->Par.me, diff);
     3748  P->R.PsiStep = 0; // make it not advance to next Psi
     3749 
     3750  //debug(P,"UpdateActualPsiNo");
     3751  UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
     3752  //debug(P,"UpdateEnergyArray");
     3753  UpdateEnergyArray(P); // shift energy values in their array by one
     3754  //debug(P,"UpdatePerturbedEnergyCalculation");
     3755  UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
     3756  EnergyAllReduce(P); // gather from all processes and sum up to total energy
     3757 
     3758  // checkout gradient
     3759  //diff = 0;
     3760  for (i=0; i< n; i+=2) {
     3761    //if ((-gradient[i/2].re != gsl_vector_get (g, i)) || (-gradient[i/2].im != gsl_vector_get (g, i+1)))  diff++;
     3762    gsl_vector_set (g, i, -gradient[i/2].re);   // real part
     3763    gsl_vector_set (g, i+1, -gradient[i/2].im); // imaginary part
     3764  }
     3765  //if (diff) fprintf(stderr,"(%i) %i differences between old and new gradient.\n", P->Par.me, diff);
     3766}
     3767
     3768/** evaluates perturbed energy functional and gradient.
     3769 * \param *x current position in functional
     3770 * \param *params void-pointer to parameter array
     3771 * \param *f pointer to energy function value on return
     3772 * \param *g array for gradient vector on return
     3773 */
     3774void perturbed_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *g) {
     3775  struct Problem *P = (struct Problem *)params;
     3776  int i, n = P->R.LevS->MaxG*2;
     3777  fftw_complex *currentPsi = P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo];
     3778  fftw_complex *gradient = P->Grad.GradientArray[ActualGradient];
     3779  //int diff = 0;
     3780  //debug(P,"fdf");
     3781  // put x into current Psi
     3782  for (i=0; i< n; i+=2) {
     3783    //if ((currentPsi[i/2].re != gsl_vector_get (x, i)) || (currentPsi[i/2].im != gsl_vector_get (x, i+1))) diff++;
     3784    currentPsi[i/2].re = gsl_vector_get (x, i);   // real part
     3785    currentPsi[i/2].im = gsl_vector_get (x, i+1); // imaginary part
     3786  }
     3787  //if (diff) fprintf(stderr,"(%i) %i differences between old and new currentPsi.\n", P->Par.me, diff);
     3788  P->R.PsiStep = 0; // make it not advance to next Psi
     3789 
     3790  //debug(P,"UpdateActualPsiNo");
     3791  UpdateActualPsiNo(P, P->R.CurrentMin); // orthogonalize
     3792  //debug(P,"UpdateEnergyArray");
     3793  UpdateEnergyArray(P); // shift energy values in their array by one
     3794  //debug(P,"UpdatePerturbedEnergyCalculation");
     3795  UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
     3796  EnergyAllReduce(P); // gather from all processes and sum up to total energy
     3797 
     3798  // checkout gradient
     3799  //diff = 0;
     3800  for (i=0; i< n; i+=2) {
     3801    //if ((-gradient[i/2].re != gsl_vector_get (g, i)) || (-gradient[i/2].im != gsl_vector_get (g, i+1)))  diff++;
     3802    gsl_vector_set (g, i, -gradient[i/2].re);   // real part
     3803    gsl_vector_set (g, i+1, -gradient[i/2].im); // imaginary part
     3804  }
     3805  //if (diff) fprintf(stderr,"(%i) %i differences between old and new gradient.\n", P->Par.me, diff);
     3806 
     3807  *f = P->Lat.E->TotalEnergy[0];   // and return evaluated functional
     3808}
     3809
     3810/* MinimisePerturbed with all the brent minimisation approach
     3811void MinimisePerturbed (struct Problem *P, int *Stop, int *SuperStop) {
     3812  struct RunStruct *R = &P->R;
     3813  struct Lattice *Lat = &P->Lat;
     3814  struct Psis *Psi = &Lat->Psi;
     3815  int type;
     3816  //int i;
     3817 
     3818  // stuff for GSL minimization
     3819  //size_t iter;
     3820  //int status, Status
     3821  int n = R->LevS->MaxG*2;
     3822  const gsl_multimin_fdfminimizer_type *T_multi;
     3823  const gsl_min_fminimizer_type *T;
     3824  gsl_multimin_fdfminimizer *s_multi;
     3825  gsl_min_fminimizer *s;
     3826  gsl_vector *x;//, *ss;
     3827  gsl_multimin_function_fdf my_func;
     3828  gsl_function F;
     3829  //fftw_complex *currentPsi;
     3830  //double a,b,m, f_m, f_a, f_b;
     3831  //double old_norm;
     3832 
     3833  my_func.f = &perturbed_f;
     3834  my_func.df = &perturbed_df;
     3835  my_func.fdf = &perturbed_fdf;
     3836  my_func.n = n;
     3837  my_func.params = P;
     3838  F.function = &perturbed_function;
     3839  F.params = P;
     3840     
     3841  x = gsl_vector_alloc (n);
     3842  //ss = gsl_vector_alloc (Psi->NoOfPsis);
     3843  T_multi = gsl_multimin_fdfminimizer_vector_bfgs;
     3844  s_multi = gsl_multimin_fdfminimizer_alloc (T_multi, n);
     3845  T = gsl_min_fminimizer_brent;
     3846  s = gsl_min_fminimizer_alloc (T);
     3847 
     3848  for (type=Perturbed_P0;type<=Perturbed_RxP2;type++) {  // go through each perturbation group separately //
     3849    *Stop=0;   // reset stop flag
     3850    fprintf(stderr,"(%i)Beginning perturbed minimisation of type %s ...\n", P->Par.me, R->MinimisationName[type]);
     3851    //OutputOrbitalPositions(P, Occupied);
     3852    R->PsiStep = R->MaxPsiStep; // reset in-Psi-minimisation-counter, so that we really advance to the next wave function
     3853    UpdateActualPsiNo(P, type); // step on to next perturbed one
     3854    fprintf(stderr, "(%i) Re-initializing perturbed psi array for type %s ", P->Par.me, R->MinimisationName[type]);
     3855    if (P->Call.ReadSrcFiles && ReadSrcPsiDensity(P,type,1, R->LevSNo)) {
     3856      SpeedMeasure(P, InitSimTime, StartTimeDo); 
     3857      fprintf(stderr,"from source file of recent calculation\n");
     3858      ReadSrcPsiDensity(P,type, 0, R->LevSNo);
     3859      ResetGramSchTagType(P, Psi, type, IsOrthogonal); // loaded values are orthonormal
     3860      SpeedMeasure(P, DensityTime, StartTimeDo);
     3861      //InitDensityCalculation(P);
     3862      SpeedMeasure(P, DensityTime, StopTimeDo);
     3863      R->OldActualLocalPsiNo = R->ActualLocalPsiNo; // needed otherwise called routines in function below crash
     3864      UpdateGramSchOldActualPsiNo(P,Psi);
     3865      InitPerturbedEnergyCalculation(P, 1);  // go through all orbitals calculate each H^{(0)}-eigenvalue, recalc HGDensity, cause InitDensityCalc zero'd it
     3866      UpdatePerturbedEnergyCalculation(P);  // H1cGradient and Gradient must be current ones
     3867      EnergyAllReduce(P);   // gather energies for minimum search
     3868      SpeedMeasure(P, InitSimTime, StopTimeDo); 
     3869    }
     3870    if (P->Call.ReadSrcFiles != 1) {
     3871      SpeedMeasure(P, InitSimTime, StartTimeDo); 
     3872      ResetGramSchTagType(P, Psi, type, NotOrthogonal); // perturbed now shall be orthonormalized
     3873      if (P->Call.ReadSrcFiles != 2) {
     3874        if (R->LevSNo == Lat->MaxLevel-1) { // is it the starting level? (see InitRunLevel())
     3875          fprintf(stderr, "randomly.\n");
     3876          InitPsisValue(P, Psi->TypeStartIndex[type], Psi->TypeStartIndex[type+1]);  // initialize perturbed array for this run
     3877        } else {
     3878          fprintf(stderr, "from source file of last level.\n");
     3879          ReadSrcPerturbedPsis(P, type);
     3880        }
     3881      }
     3882      SpeedMeasure(P, InitGramSchTime, StartTimeDo); 
     3883      GramSch(P, R->LevS, Psi, Orthogonalize);
     3884      SpeedMeasure(P, InitGramSchTime, StopTimeDo); 
     3885      SpeedMeasure(P, InitDensityTime, StartTimeDo);
     3886      //InitDensityCalculation(P);
     3887      SpeedMeasure(P, InitDensityTime, StopTimeDo);
     3888      InitPerturbedEnergyCalculation(P, 1);  // go through all orbitals calculate each H^{(0)}-eigenvalue, recalc HGDensity, cause InitDensityCalc zero'd it
     3889      R->OldActualLocalPsiNo = R->ActualLocalPsiNo; // needed otherwise called routines in function below crash
     3890      UpdateGramSchOldActualPsiNo(P,Psi);
     3891      UpdatePerturbedEnergyCalculation(P);  // H1cGradient and Gradient must be current ones
     3892      EnergyAllReduce(P);   // gather energies for minimum search
     3893      SpeedMeasure(P, InitSimTime, StopTimeDo); 
     3894      R->LevS->Step++;
     3895      EnergyOutput(P,0);
     3896      while (*Stop != 1) {
     3897        // copy current Psi into starting vector
     3898        currentPsi = R->LevS->LPsi->LocalPsi[R->ActualLocalPsiNo];
     3899        for (i=0; i< n; i+=2) {
     3900          gsl_vector_set (x, i, currentPsi[i/2].re);   // real part
     3901          gsl_vector_set (x, i+1, currentPsi[i/2].im); // imaginary part
     3902        }
     3903        gsl_multimin_fdfminimizer_set (s_multi, &my_func, x, 0.01, 1e-2); 
     3904        iter = 0;
     3905        status = 0;
     3906        do {  // look for minimum along current local psi
     3907          iter++;
     3908          status = gsl_multimin_fdfminimizer_iterate (s_multi);
     3909          MPI_Allreduce(&status, &Status, 1, MPI_INT, MPI_MAX, P->Par.comm_ST_Psi);
     3910          if (Status)
     3911            break;
     3912          status = gsl_multimin_test_gradient (s_multi->gradient, 1e-2);
     3913          MPI_Allreduce(&status, &Status, 1, MPI_INT, MPI_MAX, P->Par.comm_ST_Psi);
     3914          //if (Status == GSL_SUCCESS)
     3915            //printf ("Minimum found at:\n");
     3916          if (P->Par.me == 0) fprintf (stderr,"(%i,%i,%i)S(%i,%i,%i):\t %5d %10.5f\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, (int)iter, s_multi->f);
     3917          //TestGramSch(P,R->LevS,Psi, type); // functions are orthonormal?
     3918        } while (Status == GSL_CONTINUE && iter < 3);
     3919        // now minimize norm of currentPsi (one-dim)
     3920        if (0) {
     3921          iter = 0;
     3922        status = 0;
     3923        m = 1.;
     3924        a = MYEPSILON;
     3925        b = 100.;
     3926        f_a = perturbed_function (a, P);
     3927        f_b = perturbed_function (b, P);
     3928        f_m = perturbed_function (m, P);
     3929        //if ((f_m < f_a) && (f_m < f_b)) {
     3930          gsl_min_fminimizer_set (s, &F, m, a, b);
     3931          do {  // look for minimum along current local psi
     3932            iter++;
     3933            status = gsl_min_fminimizer_iterate (s);
     3934            m = gsl_min_fminimizer_x_minimum (s);
     3935            a = gsl_min_fminimizer_x_lower (s);
     3936            b = gsl_min_fminimizer_x_upper (s);         
     3937            status = gsl_min_test_interval (a, b, 0.001, 0.0);
     3938            if (status == GSL_SUCCESS)
     3939              printf ("Minimum found at:\n");
     3940            printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
     3941               (int) iter, a, b,
     3942               m, b - a);
     3943          } while (status == GSL_CONTINUE && iter < 100);
     3944          old_norm = GramSchGetNorm2(P,P->R.LevS,P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo]);
     3945          for (i=0; i< n; i++) {
     3946            currentPsi[i].re *= m/old_norm; // real part
     3947            currentPsi[i].im *= m/old_norm; // imaginary part
     3948          }
     3949        } else debug(P,"Norm not minimizable!");
     3950        //P->R.PsiStep = P->R.MaxPsiStep; // make it advance to next Psi
     3951        FindPerturbedMinimum(P);       
     3952        //debug(P,"UpdateActualPsiNo");
     3953        UpdateActualPsiNo(P, type); // step on to next perturbed Psi
     3954        //debug(P,"UpdateEnergyArray");
     3955        UpdateEnergyArray(P); // shift energy values in their array by one
     3956        //debug(P,"UpdatePerturbedEnergyCalculation");
     3957        UpdatePerturbedEnergyCalculation(P);  // re-calc energies (which is hopefully lower)
     3958        EnergyAllReduce(P); // gather from all processes and sum up to total energy
     3959        //ControlNativeDensity(P);  // check total density (summed up PertMixed must be zero!)
     3960        //printf ("(%i,%i,%i)S(%i,%i,%i):\t %5d %10.5f\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, (int)iter, s_multi->f);
     3961        if (*SuperStop != 1)
     3962          *SuperStop = CheckCPULIM(P);
     3963        *Stop = CalculateMinimumStop(P, *SuperStop);
     3964        P->Speed.Steps++; // step on
     3965        R->LevS->Step++;
     3966      }
     3967      // now release normalization condition and minimize wrt to norm
     3968      *Stop = 0;
     3969      while (*Stop != 1) {
     3970        currentPsi = R->LevS->LPsi->LocalPsi[R->ActualLocalPsiNo];
     3971        iter = 0;
     3972        status = 0;
     3973        m = 1.;
     3974        a = 0.001;
     3975        b = 10.;
     3976        f_a = perturbed_function (a, P);
     3977        f_b = perturbed_function (b, P);
     3978        f_m = perturbed_function (m, P);
     3979        if ((f_m < f_a) && (f_m < f_b)) {
     3980          gsl_min_fminimizer_set (s, &F, m, a, b);
     3981          do {  // look for minimum along current local psi
     3982            iter++;
     3983            status = gsl_min_fminimizer_iterate (s);
     3984            m = gsl_min_fminimizer_x_minimum (s);
     3985            a = gsl_min_fminimizer_x_lower (s);
     3986            b = gsl_min_fminimizer_x_upper (s);         
     3987            status = gsl_min_test_interval (a, b, 0.001, 0.0);
     3988            if (status == GSL_SUCCESS)
     3989              printf ("Minimum found at:\n");
     3990            printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
     3991               (int) iter, a, b,
     3992               m, b - a);
     3993          } while (status == GSL_CONTINUE && iter < 100);
     3994          old_norm = GramSchGetNorm2(P,P->R.LevS,P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo]);
     3995          for (i=0; i< n; i++) {
     3996            currentPsi[i].re *= m/old_norm; // real part
     3997            currentPsi[i].im *= m/old_norm; // imaginary part
     3998          }
     3999        }
     4000        P->R.PsiStep = P->R.MaxPsiStep; // make it advance to next Psi
     4001        //debug(P,"UpdateActualPsiNo");
     4002        UpdateActualPsiNo(P, type); // step on to next perturbed Psi
     4003        if (*SuperStop != 1)
     4004          *SuperStop = CheckCPULIM(P);
     4005        *Stop = CalculateMinimumStop(P, *SuperStop);
     4006        P->Speed.Steps++; // step on
     4007        R->LevS->Step++;
     4008      }
     4009      if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Write %s srcpsi to disk\n", P->Par.me, R->MinimisationName[type]);
     4010      OutputSrcPsiDensity(P, type);
     4011//      if (!TestReadnWriteSrcDensity(P,type))
     4012//        Error(SomeError,"TestReadnWriteSrcDensity failed!");
     4013    }
     4014
     4015    TestGramSch(P,R->LevS,Psi, type); // functions are orthonormal?
     4016    // calculate current density summands
     4017    //if (P->Call.out[StepLeaderOut]) fprintf(stderr,"(%i) Filling current density grid ...\n",P->Par.me);
     4018    SpeedMeasure(P, CurrDensTime, StartTimeDo);
     4019    if (*SuperStop != 1) {
     4020      if ((R->DoFullCurrent == 1) || ((R->DoFullCurrent == 2) && (CheckOrbitalOverlap(P) == 1))) { //test to check whether orbitals have mutual overlap and thus \\DeltaJ_{xc} must not be dropped
     4021        R->DoFullCurrent = 1; // set to 1 if it was 2 but Check...() yielded necessity
     4022        //debug(P,"Filling with Delta j ...");
     4023        //FillDeltaCurrentDensity(P);
     4024      }// else
     4025        //debug(P,"There is no overlap between orbitals.");
     4026      //debug(P,"Filling with j ...");
     4027      FillCurrentDensity(P);
     4028    }
     4029    SpeedMeasure(P, CurrDensTime, StopTimeDo);
     4030   
     4031    SetGramSchExtraPsi(P,Psi,NotUsedToOrtho); // remove extra Psis from orthogonality check
     4032    ResetGramSchTagType(P, Psi, type, NotUsedToOrtho);  // remove this group from the check for the next minimisation group as well!
     4033  }
     4034  UpdateActualPsiNo(P, Occupied); // step on back to an occupied one
     4035 
     4036  gsl_multimin_fdfminimizer_free (s_multi);
     4037  gsl_min_fminimizer_free (s);
     4038  gsl_vector_free (x);
     4039  //gsl_vector_free (ss);
     4040}
     4041*/
Note: See TracChangeset for help on using the changeset viewer.