Changeset 774ae8 for pcp

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

36f85c

Parents:

b74e5e

git-author:

Frederik Heber <heber@…> (04/18/08 15:21:30)

git-committer:

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

Message:

WriteParameters() added
OutputIonCoordinates() uses WriteParameters() and additional parameter whether its MD or StructOpt

Location:

pcp/src

Files:

• init.c (modified) (1 diff)
• init.h (modified) (1 diff)
• ions.c (modified) (1 diff)
• ions.h (modified) (1 diff)
• run.c (modified) (6 diffs)

pcp/src/init.c

@@ -1782 +1782 @@
 }
 
+/** Writes parameters as a parsable config file to \a filename.
+ * \param *P Problem at hand
+ * \param filename name of config file to be written
+ */
+void WriteParameters(struct Problem *const  P, const char *const filename) {
+  struct Lattice *Lat = &P->Lat;
+  struct RunStruct *R = &P->R;
+  struct FileData *F = &P->Files;
+  struct Ions *I = &P->Ion;
+  FILE *output;
+  double BorAng = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
+  int i, it, nr = 0;
+  if ((P->Par.me == 0) && (output = fopen(filename,"w"))) {
+    //fprintf(stderr, "(%i) Writing config file %s\n",P->Par.me, filename);
+    
+    fprintf(output,"# ParallelCarParinello - main configuration file - optimized geometry\n");
+    fprintf(output,"\n");
+    fprintf(output,"mainname\t%s\t# programm name (for runtime files)\n", F->mainname);
+    fprintf(output,"defaultpath\t%s\t# where to put files during runtime\n", F->default_path);
+    fprintf(output,"pseudopotpath\t%s\t# where to find pseudopotentials\n", F->pseudopot_path);
+    fprintf(output,"\n");
+    fprintf(output,"ProcPEGamma\t%i\t# for parallel computing: share constants\n", P->Par.proc[PEGamma]);
+    fprintf(output,"ProcPEPsi\t%i\t# for parallel computing: share wave functions\n", P->Par.proc[PEPsi]);
+    fprintf(output,"DoOutVis\t%i\t# Output data for OpenDX\n", F->DoOutVis);
+    fprintf(output,"DoOutMes\t%i\t# Output data for measurements\n", F->DoOutMes);
+    fprintf(output,"DoOutCurr\t%i\t# Ouput current density for OpenDx\n", F->DoOutCurr);
+    fprintf(output,"DoPerturbation\t%i\t# Do perturbation calculate and determine susceptibility and shielding\n", R->DoPerturbation);
+    fprintf(output,"DoFullCurrent\t%i\t# Do full perturbation\n", R->DoFullCurrent);
+    fprintf(output,"DoOutOrbitals\t%i\t# Output all orbitals on end\n", F->DoOutOrbitals);
+    fprintf(output,"CommonWannier\t%i\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center\n", R->CommonWannier);
+    fprintf(output,"SawtoothStart\t%lg\t# Absolute value for smooth transition at cell border \n", Lat->SawtoothStart);
+    fprintf(output,"VectorPlane\t%i\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot\n", R->VectorPlane);
+    fprintf(output,"VectorCut\t%lg\t# Cut plane axis value\n", R->VectorCut);
+    fprintf(output,"AddGramSch\t%i\t# Additional GramSchmidtOrtogonalization to be safe (1 - per MinStep, 2 - per PsiStep)\n", R->UseAddGramSch);
+    fprintf(output,"Seed\t\t%i\t# initial value for random seed for Psi coefficients\n", R->Seed);
+    fprintf(output,"DoBrent\t%i\t# Brent line search instead of CG with second taylor approximation\n", R->DoBrent);
+    fprintf(output,"\n");
+    fprintf(output,"MaxOuterStep\t%i\t# number of MolecularDynamics/Structure optimization steps\n", R->MaxOuterStep);
+    fprintf(output,"Deltat\t%lg\t# time in atomic time per MD step\n", R->delta_t);
+    fprintf(output,"OutVisStep\t%i\t#  Output visual data every ...th step\n", R->OutVisStep);
+    fprintf(output,"OutSrcStep\t%i\t# Output \"restart\" data every ..th step\n", R->OutSrcStep);
+    fprintf(output,"TargetTemp\t%lg\t# Target temperature in Hartree\n", R->TargetTemp);
+    fprintf(output,"Thermostat\t%s", I->ThermostatNames[I->Thermostat]);
+    switch(I->Thermostat) {
+      default:
+      case 0: // None
+        break;
+      case 1: // Woodcock
+        fprintf(output,"\t%i", R->ScaleTempStep);
+        break;
+      case 2: // Gaussian
+        fprintf(output,"\t%i", R->ScaleTempStep);
+        break;
+      case 3: // Langevin
+        fprintf(output,"\t%lg\t%lg", R->TempFrequency, R->alpha);
+        break;
+      case 4: // Berendsen
+        fprintf(output,"\t%lg", R->TempFrequency);
+        break;
+      case 5: // NoseHoover
+        fprintf(output,"\t%lg", R->HooverMass);
+        break;
+    }
+    fprintf(output,"\t# Thermostat to be used with parameters\n");
+    fprintf(output,"MaxPsiStep\t%i\t# number of Minimisation steps per state (0 - default)\n", R->MaxPsiStep);
+    fprintf(output,"EpsWannier\t%lg\t# tolerance value for spread minimisation of orbitals\n", R->EpsWannier);
+    fprintf(output,"\n");
+    fprintf(output,"# Values specifying when to stop\n");
+    fprintf(output,"MaxMinStep\t%i\t# Maximum number of steps\n", R->MaxMinStep);
+    fprintf(output,"RelEpsTotalE\t%lg\t# relative change in total energy\n", R->RelEpsTotalEnergy);
+    fprintf(output,"RelEpsKineticE\t%lg\t# relative change in kinetic energy\n", R->RelEpsKineticEnergy);
+    fprintf(output,"MaxMinStopStep\t%i\t# check every ..th steps\n", R->MaxMinStopStep);
+    fprintf(output,"MaxMinGapStopStep\t%i\t# check every ..th steps\n", R->MaxMinGapStopStep);
+    fprintf(output,"\n");
+    fprintf(output,"# Values specifying when to stop for INIT, otherwise same as above\n");
+    fprintf(output,"MaxInitMinStep\t%i\t# Maximum number of steps\n", R->MaxInitMinStep);
+    fprintf(output,"InitRelEpsTotalE\t%lg\t# relative change in total energy\n", R->InitRelEpsTotalEnergy);
+    fprintf(output,"InitRelEpsKineticE\t%lg\t# relative change in kinetic energy\n", R->InitRelEpsKineticEnergy);
+    fprintf(output,"InitMaxMinStopStep\t%i\t# check every ..th steps\n", R->InitMaxMinStopStep);
+    fprintf(output,"InitMaxMinGapStopStep\t%i\t# check every ..th steps\n", R->InitMaxMinGapStopStep);
+    fprintf(output,"\n");
+    fprintf(output,"BoxLength\t\t\t# (Length of a unit cell)\n");
+    fprintf(output,"%lg\t\n", Lat->RealBasis[0]*BorAng);
+    fprintf(output,"%lg\t%lg\t\n", Lat->RealBasis[3]*BorAng, Lat->RealBasis[4]*BorAng);
+    fprintf(output,"%lg\t%lg\t%lg\t\n", Lat->RealBasis[6]*BorAng, Lat->RealBasis[7]*BorAng, Lat->RealBasis[8]*BorAng);
+    // FIXME
+    fprintf(output,"\n");
+    fprintf(output,"ECut\t\t%lg\t# energy cutoff for discretization in Hartrees\n", Lat->ECut/(Lat->LevelSizes[0]*Lat->LevelSizes[0]));
+    fprintf(output,"MaxLevel\t%i\t# number of different levels in the code, >=2\n", P->Lat.MaxLevel);
+    fprintf(output,"Level0Factor\t%i\t# factor by which node number increases from S to 0 level\n", P->Lat.Lev0Factor);
+    fprintf(output,"RiemannTensor\t%i\t# (Use metric)\n", P->Lat.RT.Use);
+    switch (P->Lat.RT.Use) {
+      case 0: //UseNoRT
+        break;
+      case 1: // UseRT
+        fprintf(output,"RiemannLevel\t%i\t# Number of Riemann Levels\n", P->Lat.RT.RiemannLevel);
+        fprintf(output,"LevRFactor\t%i\t# factor by which node number increases from 0 to R level from\n", P->Lat.LevRFactor);
+        break;
+    }
+    fprintf(output,"PsiType\t\t%i\t# 0 - doubly occupied, 1 - SpinUp,SpinDown\n", P->Lat.Psi.Use);
+    // write out both types for easier changing afterwards
+  //  switch (PsiType) {
+  //    case 0:
+        fprintf(output,"MaxPsiDouble\t%i\t# here: specifying both maximum number of SpinUp- and -Down-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoDouble]);
+  //      break;
+  //    case 1:
+        fprintf(output,"PsiMaxNoUp\t%i\t# here: specifying maximum number of SpinUp-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoUp]);
+        fprintf(output,"PsiMaxNoDown\t%i\t# here: specifying maximum number of SpinDown-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoDown]);
+  //      break;
+  //  }
+    fprintf(output,"AddPsis\t\t%i\t# Additional unoccupied Psis for bandgap determination\n", P->Lat.Psi.GlobalNo[PsiMaxAdd]);
+    fprintf(output,"\n");
+    fprintf(output,"RCut\t\t%lg\t# R-cut for the ewald summation\n", I->R_cut);
+    fprintf(output,"StructOpt\t%i\t# Do structure optimization beforehand\n", I->StructOpt);
+    fprintf(output,"IsAngstroem\t%i\t# 0 - Bohr, 1 - Angstroem\n", I->IsAngstroem);
+    fprintf(output,"RelativeCoord\t0\t# whether ion coordinates are relative (1) or absolute (0)\n");
+    fprintf(output,"MaxTypes\t%i\t# maximum number of different ion types\n", I->Max_Types);
+    fprintf(output,"\n");
+    // now elements ...
+    fprintf(output,"# Ion type data (PP = PseudoPotential, Z = atomic number)\n");
+    fprintf(output,"#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\tName\tSymbol\n");
+    for (i=0; i < I->Max_Types; i++)
+      fprintf(output,"Ion_Type%i\t%i\t%i\t%lg\t%i\t%i\t%lg\t%s\t%s\n", i+1, I->I[i].Max_IonsOfType, I->I[i].Z, I->I[i].rgauss, I->I[i].l_max, I->I[i].l_loc, I->I[i].IonMass/Units2Electronmass, &I->I[i].Name[0], &I->I[i].Symbol[0]);
+    // ... and each ion coordination
+    fprintf(output,"#Ion_TypeNr._Nr.R[0]\tR[1]\tR[2]\tMoveType (0 MoveIon, 1 FixedIon) U[0]\tU[1]\tU[2] (velocities in %s)\n", (I->IsAngstroem) ? "Angstroem" : "atomic units");
+    for (i=0; i < I->Max_Types; i++)
+      for (it=0;it<I->I[i].Max_IonsOfType;it++) {
+        fprintf(output,"Ion_Type%i_%i\t%2.9f\t%2.9f\t%2.9f\t%i\t%2.9f\t%2.9f\t%2.9f\t# Number in molecule %i\n", i+1, it+1, I->I[i].R[0+NDIM*it]*BorAng, I->I[i].R[1+NDIM*it]*BorAng, I->I[i].R[2+NDIM*it]*BorAng, I->I[i].IMT[it], I->I[i].U[0+NDIM*it]*BorAng, I->I[i].U[1+NDIM*it]*BorAng, I->I[i].U[2+NDIM*it]*BorAng, ++nr);
+      }
+    fflush(output);
+    fclose(output);
+  }
+}
+
 /** Main Initialization.
  * Massive calling of initializing super-functions init...():

pcp/src/init.h

@@ -31 +31 @@
 int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, int repetition, enum necessity critical);
 void ReadParameters(struct Problem *const  P, const char *const filename);
+void WriteParameters(struct Problem *const  P, const char *const filename);
 
 void Init(struct Problem *const P);

pcp/src/ions.c

@@ -959 +959 @@
 /** Print coordinates of all ions to stdout.
  * \param *P Problem at hand
- */
-void OutputIonCoordinates(struct Problem *P)
-{
-  struct Ions *I = &P->Ion;
-  int is, ia;
-  if (P->Par.me == 0) { 
-    fprintf(stderr, "(%i) ======= Updated Ion Coordinates =========\n",P->Par.me);
+ * \param actual (1 - don't at current RunStruct#OuterStep, 0 - do)
+ */
+void OutputIonCoordinates(struct Problem *P, int actual)
+{
+  //struct RunStruct *R = &P->R;
+  struct Ions *I = &P->Ion;
+  char filename[255];
+  FILE *output;
+  int is, ia, nr = 0;
+  double Bohr = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
+/*  if (P->Par.me == 0 && P->Call.out[ReadOut]) { 
+    fprintf(stderr, "(%i) ======= Differential Ion Coordinates =========\n",P->Par.me);
     for (is=0; is < I->Max_Types; is++)
       for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
         //fprintf(stderr, "(%i) R[%i/%i][%i/%i] = (%e,%e,%e)\n", P->Par.me, is, I->Max_Types, ia, I->I[is].Max_IonsOfType, I->I[is].R[NDIM*ia+0],I->I[is].R[NDIM*ia+1],I->I[is].R[NDIM*ia+2]);  
-        fprintf(stderr, "Ion_Type%i_%i\t%.6f\t%.6f\t%.6f\t0\t# Atom from StructOpt\n", is+1, ia+1, I->I[is].R[NDIM*ia+0],I->I[is].R[NDIM*ia+1],I->I[is].R[NDIM*ia+2]);
+        fprintf(stderr, "Ion_Type%i_%i\t%.6f\t%.6f\t%.6f\t0\t# Atom from StructOpt\n", is+1, ia+1, I->I[is].R[NDIM*ia+0]-I->I[is].R_old[NDIM*ia+0],I->I[is].R[NDIM*ia+1]-I->I[is].R_old[NDIM*ia+1],I->I[is].R[NDIM*ia+2]-I->I[is].R_old[NDIM*ia+2]);
     fprintf(stderr, "(%i) =========================================\n",P->Par.me);
+  }*/
+  if (actual)
+    sprintf(filename, "%s.MD", P->Call.MainParameterFile); 
+  else
+    sprintf(filename, "%s.opt", P->Call.MainParameterFile);
+  if (((actual) && (P->R.OuterStep <= 0)) || ((!actual) && (P->R.StructOptStep == 1))) { // on first step write complete config file 
+    WriteParameters(P, filename);
+  } else { // afterwards simply add the current ion coordinates as constrained steps
+    if ((P->Par.me == 0) && (output = fopen(filename,"a"))) {
+      fprintf(output,"# ====== MD step %d ========= \n", (actual) ? P->R.OuterStep : P->R.StructOptStep);
+      for (is=0; is < I->Max_Types; is++)
+        for (ia=0;ia<I->I[is].Max_IonsOfType;ia++) {
+          fprintf(output,"Ion_Type%i_%i\t%2.9f\t%2.9f\t%2.9f\t%i\t%2.9f\t%2.9f\t%2.9f\t# Number in molecule %i\n", is+1, ia+1, I->I[is].R[0+NDIM*ia]*Bohr, I->I[is].R[1+NDIM*ia]*Bohr, I->I[is].R[2+NDIM*ia]*Bohr, I->I[is].IMT[ia], I->I[is].U[0+NDIM*ia]*Bohr, I->I[is].U[1+NDIM*ia]*Bohr, I->I[is].U[2+NDIM*ia]*Bohr, ++nr);
+        }
+      fflush(output);
+    }
   }
 }

pcp/src/ions.h

@@ -128 +128 @@
 void OutputIonForce(struct Problem *P);
 void ParseIonForce(struct Problem *P);
-void OutputIonCoordinates(struct Problem *P);
+void OutputIonCoordinates(struct Problem *P, int actual);
 void UpdateIons(struct Problem *P);
 void UpdateIonsR(struct Problem *P);

pcp/src/run.c

@@ -1429 +1429 @@
     //gsl_vector_fprintf(stderr, s->dx, "%lg");
     OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
-    //OutputIonCoordinates(P, 0);
+    OutputIonCoordinates(P, 0);
     P->R.StructOptStep++;
   } while ((status == GSL_CONTINUE) && (P->R.StructOptStep < P->R.MaxStructOptStep));
@@ -1499 +1499 @@
     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);
+    OutputIonCoordinates(P, 0);
     P->R.StructOptStep++;
   } while ((status == GSL_CONTINUE) && (Run->OuterStep < Run->MaxOuterStep));
@@ -1742 +1742 @@
   SpeedMeasure(P, InitSimTime, StopTimeDo);
 
-  //OutputIonCoordinates(P, 1);
+  OutputIonCoordinates(P, 1);
   OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
   OutputIonForce(P);
@@ -1755 +1755 @@
     while ((R->MeanForce[0] > 1e-4) && (R->StructOptStep < R->MaxStructOptStep)) {
       R->StructOptStep++; 
-      //OutputIonCoordinates(P, 1);
+      OutputIonCoordinates(P, 1);
       UpdateIons(P); 
       UpdateWaveAfterIonMove(P);
@@ -1773 +1773 @@
       if (P->Par.me == 0) fprintf(stderr,"(%i) Mean force is %lg\n", P->Par.me, R->MeanForce[0]); 
     }
-    //OutputIonCoordinates(P, 1);
+    OutputIonCoordinates(P, 1);
   }
   if (I->StructOpt && !OuterStop) {
@@ -1799 +1799 @@
 
     UpdateIonsR(P);
-    //OutputIonCoordinates(P, 1);
+    OutputIonCoordinates(P, 1);
 
     UpdateWaveAfterIonMove(P);