Changeset 2026d4 for pcp

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

27a5bf

Parents:

963310a

git-author:

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

git-committer:

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

Message:

StructOpt_func() and CheckForChangedPositions() added

File:

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

pcp/src/run.c

@@ -1153 +1153 @@
           SpeedMeasure(P, SimTime, StopTimeDo);
           ControlNativeDensity(P);
-          TestGramSch(P,LevS,Psi, Occupied);   
-    // necessary for correct ionic forces ...
-    SpeedMeasure(P, LocFTime, StartTimeDo);
-    CalculateIonLocalForce(P);
-    SpeedMeasure(P, LocFTime, StopTimeDo);
-    SpeedMeasure(P, NonLocFTime, StartTimeDo);
-    CalculateIonNonLocalForce(P);
-    SpeedMeasure(P, NonLocFTime, StopTimeDo);
-    CalculateEwald(P, 1);
-    CalculateIonForce(P);
+          TestGramSch(P,LevS,Psi, Occupied);
+          // necessary for correct ionic forces ...
+                SpeedMeasure(P, LocFTime, StartTimeDo);
+                CalculateIonLocalForce(P);
+                SpeedMeasure(P, LocFTime, StopTimeDo);
+                SpeedMeasure(P, NonLocFTime, StartTimeDo);
+                CalculateIonNonLocalForce(P);
+                SpeedMeasure(P, NonLocFTime, StopTimeDo);
+                CalculateEwald(P, 1);
+                CalculateIonForce(P);
   }
   if (P->Call.ForcesFile != NULL) { // if we parse forces from file, values are written over (in case of DoOutVis)
@@ -1168 +1168 @@
     //ParseIonForce(P); 
     //CalculateIonForce(P);
-  } 
+  }
   CorrectForces(P);
   // ... on output of densities
@@ -1177 +1177 @@
 
   //OutputNorm(P);
-  //fprintf(stderr,"(%i) DoubleG: %p, CArray[22]: %p, OldLocalPsi: %p\n", P->Par.me, R->LevS->DoubleG, R->Lev0->Dens->DensityCArray[22], R->LevS->LPsi->OldLocalPsi); 
+  //fprintf(stderr,"(%i) DoubleG: %p, CArray[22]: %p, OldLocalPsi: %p\n", P->Par.me, R->LevS->DoubleG, R->Lev0->Dens->DensityCArray[22], R->LevS->LPsi->OldLocalPsi);
   //OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
-  /*TestGramSch(P,LevS,Psi, &P->Lat.Psi); */
+  /*TestGramSch(P, R->LevS, &P->Lat.Psi); */
   GetOuterStop(P);
   //fprintf(stderr,"(%i) DoubleG: %p, CArray[22]: %p, OldLocalPsi: %p\n", P->Par.me, R->LevS->DoubleG, R->Lev0->Dens->DensityCArray[22], R->LevS->LPsi->OldLocalPsi); 
   if (SuperStop) OuterStop = 1;
   return OuterStop;
+}
+
+/** Checks whether the given positions \a *v have changed wrt stored in IonData structure.
+ * \param *P Problem at hand
+ * \param *v gsl_vector storing new positions
+ */
+int CheckForChangedPositions(struct Problem *P, const gsl_vector *v)
+{
+  struct Ions *I = &P->Ion;
+  int is,ia,k, index=0;
+  int diff = 0;
+  double *R_ion;
+  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];
+      for (k=0;k<NDIM;k++) { // for all dimensions
+        if (fabs(R_ion[k] - gsl_vector_get (v, index++)) > MYEPSILON)
+          diff++;
+      }      
+    }
+  return diff;  
+}
+
+/** Wrapper for CalculateForce() for simplex minimisation of total energy.
+ * \param *v vector with degrees of freedom
+ * \param *params additional arguments, here Problem at hand
+ */
+double StructOpt_func(const gsl_vector *v, void *params)
+{
+  struct Problem *P = (struct Problem *)params;
+  struct RunStruct *R = &P->R;
+  struct Ions *I = &P->Ion;
+  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 = 0;
+  int OuterStop;
+  double diff = 0., tmp;
+  debug (P, "StructOpt_func");
+  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);*/
+    }
+  }
+  if (P->Par.me == 0) fprintf(stderr,"(%i) TE %e\n",P->Par.me, E->TotalEnergy[0]);
+  return E->TotalEnergy[0];
 }