Changeset 27a5bf

Timestamp:

Apr 21, 2008, 2:19:23 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

9b2f5d

Parents:

2026d4

git-author:

Frederik Heber <heber@…> (04/18/08 14:15:13)

git-committer:

Frederik Heber <heber@…> (04/21/08 14:19:23)

Message:

-StructOpt_...(): Huge changes (from newest svn version)
-CalculateMD(): Huge changes to get MD smooth and going

File:

• pcp/src/run.c (modified) (11 diffs)

pcp/src/run.c ¶

@@ -1281 +1281 @@
   struct Energy *E = P->Lat.E;
   int i;
-  double *R_ion, *R_old, *R_old_old, *FIon;
-  double norm = 0.;
-  int is,ia,k,index;
+  double *R_ion, *R_old, *R_old_old;//, *FIon;
+  //double norm = 0.;
+  int is,ia,k,index = 0;
   int OuterStop;
-  // update ion positions from vector coordinates
-  index=0;
-  for (is=0; is < I->Max_Types; is++) // for all elements
-    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
-      R_ion = &I->I[is].R[NDIM*ia];
-      R_old = &I->I[is].R_old[NDIM*ia];
-      R_old_old = &I->I[is].R_old_old[NDIM*ia];
-      for (k=0;k<NDIM;k++) { // for all dimensions
-        R_old_old[k] = R_old[k];
-        R_old[k] = R_ion[k];
-        R_ion[k] = gsl_vector_get (v, index++);        
-      }      
-    }  
-  // recalculate ionic forces (do electronic minimisation)
-  R->OuterStep++;
-  if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Commencing Fletcher-Reeves step %i ... \n",P->Par.me, R->OuterStep);
-  R->NewRStep++;
-  OutputIonCoordinates(P);
-  UpdateWaveAfterIonMove(P);
-  for (i=MAXOLD-1; i > 0; i--)    // store away old energies
-    E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
-  UpdateToNewWaves(P);
-  E->TotalEnergyOuter[0] = E->TotalEnergy[0];
-  OuterStop = CalculateForce(P);
-  UpdateIonsU(P);
-  CorrectVelocity(P);
-  CalculateEnergyIonsU(P);
-/*  if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
-    ScaleTemp(P);*/
-  if ((R->OuterStep-1) % P->R.OutSrcStep == 0) 
-    OutputVisSrcFiles(P, Occupied);
-  if ((R->OuterStep-1) % P->R.OutVisStep == 0) {      
-/*      // recalculate density for the specific wave function ...
-      CalculateOneDensityR(Lat, LevS, Dens0, PsiDat, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
-      // ... and output (wherein ActualDensity is used instead of TotalDensity)
-      OutputVis(P);
-      OutputIonForce(P);
-      EnergyOutput(P, 1);*/
-  }
-  // sum up mean force
-  for (is=0; is < I->Max_Types; is++)
-    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
-      FIon = &I->I[is].FIon[NDIM*ia];
-      norm += sqrt(RSP3(FIon,FIon));
-    }
-  if (P->Par.me == 0) fprintf(stderr,"(%i) Mean Force over all Ions %e\n",P->Par.me, norm);
-  return norm;
+  double diff = 0., tmp;
+  //debug (P, "StructOpt_f");
+  if (CheckForChangedPositions(P,v)) {
+    // update ion positions from vector coordinates
+    for (is=0; is < I->Max_Types; is++) // for all elements
+      for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
+        R_ion = &I->I[is].R[NDIM*ia];
+        R_old = &I->I[is].R_old[NDIM*ia];
+        R_old_old = &I->I[is].R_old_old[NDIM*ia];
+        tmp = 0.;
+        for (k=0;k<NDIM;k++) { // for all dimensions
+          R_old_old[k] = R_old[k];
+          R_old[k] = R_ion[k];
+          tmp += (R_ion[k]-gsl_vector_get (v, index))*(R_ion[k]-gsl_vector_get (v, index));
+          R_ion[k] = gsl_vector_get (v, index++);        
+        }
+        diff += sqrt(tmp);
+      }  
+    if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Summed Difference to former position %lg\n", P->Par.me, diff);
+    // recalculate ionic forces (do electronic minimisation)
+    //R->OuterStep++;
+    R->NewRStep++;
+    UpdateWaveAfterIonMove(P);
+    for (i=MAXOLD-1; i > 0; i--)    // store away old energies
+      E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
+    UpdateToNewWaves(P);
+    E->TotalEnergyOuter[0] = E->TotalEnergy[0];
+    OuterStop = CalculateForce(P);
+    //UpdateIonsU(P);
+    //CorrectVelocity(P);
+    //CalculateEnergyIonsU(P);
+  /*  if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
+      ScaleTemp(P);*/
+    if ((R->OuterStep-1) % P->R.OutSrcStep == 0) 
+      OutputVisSrcFiles(P, Occupied);
+    /*if ((R->OuterStep-1) % P->R.OutVisStep == 0) {      
+        // recalculate density for the specific wave function ...
+        CalculateOneDensityR(Lat, LevS, Dens0, PsiDat, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
+        // ... and output (wherein ActualDensity is used instead of TotalDensity)
+        OutputVis(P);
+        OutputIonForce(P);
+        EnergyOutput(P, 1);
+    }*/
+  }
+  GetOuterStop(P);
+  //if (P->Call.out[LeaderOut] && (P->Par.me == 0)) fprintf(stderr,"(%i) Absolute Force summed over all Ions %e\n",P->Par.me, norm);
+  return R->MeanForce[0];
+  //if (P->Call.out[LeaderOut] && (P->Par.me == 0)) fprintf(stderr,"(%i) Struct_optf returning: %lg\n",P->Par.me,E->TotalEnergy[0]); 
+  //return E->TotalEnergy[0];
 }
 
@@ -1339 +1341 @@
   struct Ions *I = &P->Ion;
   double *FIon;
-  int is,ia,k, index=0;  
+  int is,ia,k, index=0;
+  //debug (P, "StructOpt_df");
+  // look through coordinate vector if positions have changed sind last StructOpt_f call
+  if (CheckForChangedPositions(P,v)) {// if so, recalc to update forces
+    debug (P, "Calling StructOpt_f to update");
+    StructOpt_f(v, params);
+  }
   for (is=0; is < I->Max_Types; is++) // for all elements
     for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
@@ -1346 +1354 @@
         gsl_vector_set (df, index++, FIon[k]);
       }      
-    }   
+    } 
+  if (P->Call.out[LeaderOut] && (P->Par.me == 0)) {
+    fprintf(stderr,"(%i) Struct_Optdf returning",P->Par.me); 
+    gsl_vector_fprintf(stderr, df, "%lg");
+  } 
 }
 
@@ -1362 +1374 @@
 void UpdateIon_PRCG(struct Problem *P)
 {
-  struct RunStruct *Run = &P->R;
+  //struct RunStruct *Run = &P->R;
   struct Ions *I = &P->Ion;
   size_t np = NDIM*I->Max_TotalIons;  // d.o.f = number of ions times number of dimensions
@@ -1378 +1390 @@
   /* Starting point */
   x = gsl_vector_alloc (np);
+  //fprintf(stderr,"(%i) d.o.f. = %i\n", P->Par.me, (int)np);
 
   index=0;
@@ -1397 +1410 @@
   s = gsl_multimin_fdfminimizer_alloc (T, np);
 
-  gsl_multimin_fdfminimizer_set (s, &minex_func, x, 0.1, 1e-4);
-
+  gsl_multimin_fdfminimizer_set (s, &minex_func, x, 0.1, 0.001);
+
+  fprintf(stderr,"(%i) Commencing Structure optimization with PRCG: dof %d\n", P->Par.me,(int)np);
   do {
     iter++;
@@ -1406 +1420 @@
       break;
     
-    status = gsl_multimin_test_gradient (s->gradient, 1e-4);
+    status = gsl_multimin_test_gradient (s->gradient, 1e-2);
        
     if (status == GSL_SUCCESS)
@@ -1413 +1427 @@
     //if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Commencing '%s' step %i ... \n",P->Par.me, gsl_multimin_fdfminimizer_name(s), P->R.StructOptStep);
     if ((P->Call.out[NormalOut]) && (P->Par.me == 0)) fprintf (stderr, "(%i) %5d %10.5f\n", P->Par.me, (int)iter, s->f);
-  } while (status == GSL_CONTINUE && iter < Run->MaxOuterStep);
+    //gsl_vector_fprintf(stderr, s->dx, "%lg");
+    OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
+    //OutputIonCoordinates(P, 0);
+    P->R.StructOptStep++;
+  } while ((status == GSL_CONTINUE) && (P->R.StructOptStep < P->R.MaxStructOptStep));
 
   gsl_vector_free(x);
   gsl_multimin_fdfminimizer_free (s);
+}
+
+/** Simplex implementation for the structure optimization.
+ * We follow the example from the GSL manual.
+ * \param *P Problem at hand
+ */
+void UpdateIon_Simplex(struct Problem *P)
+{
+  struct RunStruct *Run = &P->R;
+  struct Ions *I = &P->Ion;
+  size_t np = NDIM*I->Max_TotalIons;  // d.o.f = number of ions times number of dimensions
+  int is, ia, k, index;
+  double *R;
+ 
+  const gsl_multimin_fminimizer_type *T;
+  gsl_multimin_fminimizer *s;
+  gsl_vector *x, *ss;
+  gsl_multimin_function minex_func;
+
+  size_t iter = 0;
+  int status;
+  double size;
+
+  ss = gsl_vector_alloc (np);
+  gsl_vector_set_all(ss, .2);
+  /* Starting point */
+  x = gsl_vector_alloc (np);
+  //fprintf(stderr,"(%i) d.o.f. = %i\n", P->Par.me, (int)np);
+
+  index=0;
+  for (is=0; is < I->Max_Types; is++) // for all elements
+    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {  // for all ions of element
+      R = &I->I[is].R[NDIM*ia];
+      for (k=0;k<NDIM;k++) // for all dimensions
+        gsl_vector_set (x, index++, R[k]);      
+    }
+
+  /* Initialize method and iterate */
+  minex_func.f = &StructOpt_f;
+  minex_func.n = np;
+  minex_func.params = (void *)P;
+  
+  T = gsl_multimin_fminimizer_nmsimplex;
+  s = gsl_multimin_fminimizer_alloc (T, np);
+
+  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
+
+  fprintf(stderr,"(%i) Commencing Structure optimization with NM simplex: dof %d\n", P->Par.me, (int)np);
+  do {
+    iter++;
+    status = gsl_multimin_fminimizer_iterate(s);
+
+    if (status)
+      break;
+    
+    size = gsl_multimin_fminimizer_size (s);
+    status = gsl_multimin_test_size (size, 1e-4);
+       
+    if (status == GSL_SUCCESS)
+      if (P->Par.me == 0) fprintf (stderr,"(%i) converged to minimum at\n", P->Par.me);
+
+    if (P->Call.out[MinOut]) fprintf(stderr,"(%i) Commencing '%s' step %i ... \n",P->Par.me, gsl_multimin_fminimizer_name(s), P->R.StructOptStep);
+    if ((P->Call.out[MinOut]) && (P->Par.me == 0)) fprintf (stderr, "(%i) %5d %10.5f %10.5f\n", P->Par.me, (int)iter, s->fval, size);
+    OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
+    //OutputIonCoordinates(P, 0);
+    P->R.StructOptStep++;
+  } while ((status == GSL_CONTINUE) && (Run->OuterStep < Run->MaxOuterStep));
+
+  gsl_vector_free(x);
+  gsl_vector_free(ss);
+  gsl_multimin_fminimizer_free (s);
 }
 
@@ -1638 +1727 @@
   struct RunStruct *R = &P->R;
   struct Ions *I = &P->Ion;
+  struct Energy *E = P->Lat.E;
   int OuterStop = 0;
+  int i;
+  
   SpeedMeasure(P, SimTime, StartTimeDo);
+  // initial calculations (bring density on BO surface and output start energies, coordinates, densities, ...)
   SpeedMeasure(P, InitSimTime, StartTimeDo);
   R->OuterStep = 0;
@@ -1645 +1738 @@
   CalculateEnergyIonsU(P);
   OuterStop = CalculateForce(P);
-  R->OuterStep++;
+  //R->OuterStep++;
   P->Speed.InitSteps++;
   SpeedMeasure(P, InitSimTime, StopTimeDo);
+
+  //OutputIonCoordinates(P, 1);
+  OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
   OutputIonForce(P);
   EnergyOutput(P, 1);
-  if (R->MaxOuterStep > 0) {
-    debug(P,"Commencing Fletcher-Reeves minimisation on ionic structure ...");
-    UpdateIon_PRCG(P);
+  
+  // if desired perform beforehand a structure relaxation/optimization
+  if (I->StructOpt) {
+    debug(P,"Commencing minimisation on ionic structure ...");
+    R->StructOptStep = 0;
+    //UpdateIon_PRCG(P);
+    //UpdateIon_Simplex(P);
+    while ((R->MeanForce[0] > 1e-4) && (R->StructOptStep < R->MaxStructOptStep)) {
+      R->StructOptStep++; 
+      //OutputIonCoordinates(P, 1);
+      UpdateIons(P); 
+      UpdateWaveAfterIonMove(P);
+      for (i=MAXOLD-1; i > 0; i--)    // store away old energies
+        E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
+      UpdateToNewWaves(P);
+      E->TotalEnergyOuter[0] = E->TotalEnergy[0]; 
+      OuterStop = CalculateForce(P);
+      CalculateEnergyIonsU(P);
+      if ((R->StructOptStep-1) % P->R.OutSrcStep == 0) 
+        OutputVisSrcFiles(P, Occupied);
+      if ((R->StructOptStep-1) % P->R.OutVisStep == 0) {      
+        OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
+        OutputIonForce(P);
+        EnergyOutput(P, 1);
+      }
+      if (P->Par.me == 0) fprintf(stderr,"(%i) Mean force is %lg\n", P->Par.me, R->MeanForce[0]); 
+    }
+    //OutputIonCoordinates(P, 1);
   }
   if (I->StructOpt && !OuterStop) {
@@ -1658 +1779 @@
     OuterStop = CalculateForce(P);
   }
- /*  while (!OuterStop && R->OuterStep <= R->MaxOuterStep) {
+  
+  // and now begin with the molecular dynamics simulation
+  debug(P,"Commencing MD simulation ...");
+  while (!OuterStop && R->OuterStep < R->MaxOuterStep) {
     R->OuterStep++;
-    if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Commencing MD steps %i ... \n",P->Par.me, R->OuterStep);
-    P->R.t += P->R.delta_t;     // increase current time by delta_t
+    if (P->Par.me == 0) {
+      if (R->OuterStep > 1) fprintf(stderr,"\b\b\b\b\b\b\b\b\b\b\b\b");
+      fprintf(stderr,"Time: %f fs\r", R->t*Atomictime2Femtoseconds);
+      fflush(stderr);
+    }
+    OuterStop = CalculateForce(P);
+    P->R.t += P->R.delta_t;   // increase current time by delta_t
     R->NewRStep++;
-    if (P->Ion.StructOpt == 1) {
-      UpdateIons(P);
-      OutputIonCoordinates(P);
-    } else { 
-      UpdateIonsR(P);
-    }
+
+    UpdateIonsU(P);
+    CorrectVelocity(P);
+    Thermostats(P, I->Thermostat);
+    CalculateEnergyIonsU(P);
+
+    UpdateIonsR(P);
+    //OutputIonCoordinates(P, 1);
+
     UpdateWaveAfterIonMove(P);
-    for (i=MAXOLD-1; i > 0; i--)                // store away old energies
+    for (i=MAXOLD-1; i > 0; i--)    // store away old energies
       E->TotalEnergyOuter[i] = E->TotalEnergyOuter[i-1];  
     UpdateToNewWaves(P);
     E->TotalEnergyOuter[0] = E->TotalEnergy[0]; 
-    OuterStop = CalculateForce(P);
-    UpdateIonsU(P);
-    CorrectVelocity(P);
-    CalculateEnergyIonsU(P);
-    if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
-      ScaleTemp(P);
+    //if ((P->R.ScaleTempStep > 0) && ((R->OuterStep-1) % P->R.ScaleTempStep == 0))
+    //  ScaleTemp(P);
     if ((R->OuterStep-1) % P->R.OutSrcStep == 0) 
       OutputVisSrcFiles(P, Occupied);
     if ((R->OuterStep-1) % P->R.OutVisStep == 0) {      
-      // recalculate density for the specific wave function ...
-      //CalculateOneDensityR(Lat, LevS, Dens0, PsiDat, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
-      // ... and output (wherein ActualDensity is used instead of TotalDensity)
-      //OutputVis(P);
-      //OutputIonForce(P);
-      //EnergyOutput(P, 1);
-    }
-  }*/
+      OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
+      OutputIonForce(P);
+      EnergyOutput(P, 1);
+    }
+    ResetForces(P);
+  }
   SpeedMeasure(P, SimTime, StopTimeDo);
-  // hack to output each orbital before and after spread-minimisation
-/*  if (P->Files.MeOutVis) P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+P->Lat.Psi.MaxPsiOfType*2),"OutputVis");
-  OutputVisAllOrbital(P, 0, 2, Occupied);
-  CalculateHamiltonian(P);
-  if (P->Files.MeOutVis) P->Files.OutVisStep -= (P->Lat.Psi.MaxPsiOfType)*2;
-  OutputVisAllOrbital(P, 1, 2, Occupied);*/  
   CloseOutputFiles(P);
 }