Changeset b74e5e for pcp/src

Timestamp:

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

Author:

Frederik Heber <heber@…>

Children:

774ae8

Parents:

961b75

git-author:

Frederik Heber <heber@…> (04/18/08 15:17:03)

git-committer:

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

Message:

InitThermostats() added and included in IonsInitRead()

Location:

pcp/src

Files:

• ions.c (modified) (3 diffs)
• ions.h (modified) (1 diff)

pcp/src/ions.c

@@ -50 +50 @@
 #include "ions.h"
 #include "init.h"
+
+/** Readin of Thermostat related values from parameter file.
+ * \param *P Problem at hand
+ * \param *source parameter file
+ */
+void InitThermostats(struct Problem *P, FILE *source)
+{
+  struct FileData *F = &P->Files;
+  struct Ions *I = &P->Ion;
+  char *thermo = MallocString(12, "IonsInitRead: thermo");
+  int j;
+
+  // initialize the naming stuff and all
+  I->ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
+  I->ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
+  for (j=0;j<MaxThermostats;j++)
+    I->ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
+  strcpy(I->ThermostatNames[0],"None");
+  I->ThermostatImplemented[0] = 1;
+  strcpy(I->ThermostatNames[1],"Woodcock"); 
+  I->ThermostatImplemented[1] = 1;
+  strcpy(I->ThermostatNames[2],"Gaussian"); 
+  I->ThermostatImplemented[2] = 1;
+  strcpy(I->ThermostatNames[3],"Langevin"); 
+  I->ThermostatImplemented[3] = 1;
+  strcpy(I->ThermostatNames[4],"Berendsen"); 
+  I->ThermostatImplemented[4] = 1;
+  strcpy(I->ThermostatNames[5],"NoseHoover"); 
+  I->ThermostatImplemented[5] = 1;
+  I->Thermostat = -1;
+  if (F->MeOutMes) fprintf(F->TemperatureFile, "# %s\t%s\t%s --- ","Step", "ActualTemp(1)", "ActualTemp(2)");
+
+  // read desired Thermostat from file along with needed additional parameters
+  if (ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
+    if (strcmp(thermo, I->ThermostatNames[0]) == 0) { // None
+      if (I->ThermostatImplemented[0] == 1) {
+        I->Thermostat = None;
+      } else {
+        if (P->Par.me == 0)
+          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[0]);
+        I->Thermostat = None;
+      }
+    } else if (strcmp(thermo, I->ThermostatNames[1]) == 0) { // Woodcock
+      if (I->ThermostatImplemented[1] == 1) {
+        I->Thermostat = Woodcock;
+        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, int_type, &P->R.ScaleTempStep, 1, critical); // read scaling frequency
+      } else {
+        if (P->Par.me == 0)
+          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[1]);
+        I->Thermostat = None;
+      }
+    } else if (strcmp(thermo, I->ThermostatNames[2]) == 0) { // Gaussian
+      if (I->ThermostatImplemented[2] == 1) {
+        I->Thermostat = Gaussian;
+        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, int_type, &P->R.ScaleTempStep, 1, critical); // read collision rate
+      } else {
+        if (P->Par.me == 0)
+          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[2]);
+        I->Thermostat = None;
+      }
+    } else if (strcmp(thermo, I->ThermostatNames[3]) == 0) { // Langevin
+      if (I->ThermostatImplemented[3] == 1) {
+        I->Thermostat = Langevin;
+        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.TempFrequency, 1, critical); // read gamma
+        if (ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 3, 1, double_type, &P->R.alpha, 1, optional)) {
+          if (P->Par.me == 0) 
+            fprintf(stderr,"(%i) Extended Stochastic Thermostat detected with interpolation coefficient %lg\n", P->Par.me, P->R.alpha);
+        } else {
+          P->R.alpha = 1.;
+        }
+      } else {
+        if (P->Par.me == 0)
+          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[3]);
+        I->Thermostat = None;
+      }
+    } else if (strcmp(thermo, I->ThermostatNames[4]) == 0) { // Berendsen
+      if (I->ThermostatImplemented[4] == 1) {
+        I->Thermostat = Berendsen;
+        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.TempFrequency, 1, critical); // read \tau_T
+      } else {
+        if (P->Par.me == 0)
+          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[4]);
+        I->Thermostat = None;
+      }
+    } else if (strcmp(thermo, I->ThermostatNames[5]) == 0) { // Nose-Hoover
+      if (I->ThermostatImplemented[5] == 1) {
+        I->Thermostat = NoseHoover;
+        ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.HooverMass, 1, critical); // read Hoovermass
+        P->R.alpha = 0.;
+      } else {
+        if (P->Par.me == 0)
+          fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[5]);
+        I->Thermostat = None;
+      }
+    }
+    if (I->Thermostat == -1) {
+      if (P->Par.me == 0) 
+        fprintf(stderr,"(%i) Warning: thermostat name was not understood!\n", P->Par.me);
+      I->Thermostat = None;
+    }
+  } else {
+    if ((P->R.MaxOuterStep > 0) && (P->R.TargetTemp != 0)) 
+      debug(P, "No thermostat chosen despite finite temperature MD, falling back to None.");
+    I->Thermostat = None;
+  }
+  if (F->MeOutMes) fprintf(F->TemperatureFile, "%s Thermostat\n",I->ThermostatNames[(int)I->Thermostat]); 
+  Free(thermo, "InitThermostats: thermo");
+}
 
 /** Parses for a given keyword the ion coordinates and velocites.
@@ -179 +287 @@
     I->StructOpt = 0;
   
-  //InitThermostats(P, source);
+  InitThermostats(P, source);
   
   /* Ions Data */
@@ -515 +623 @@
   fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
 } 
+
+/** Reads Forces from CallOptions::ForcesFile.
+ * goes through all Ions per IonType and each dimension of #NDIM and parses in one line:
+ * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL, 
+ * non-local IonType::FIonNL, ewald force IonType::FEwald
+ * \param *P Problem at hand
+ */
+void ParseIonForce(struct Problem *P)
+{
+        //struct RunStruct *Run = &P->R;
+  struct Ions *I = &P->Ion;
+  FILE *finput = NULL;
+        char line[1024];
+  int i,j;
+  double i1, j1;
+  double R[NDIM];
+  
+  if (P->Call.ForcesFile == NULL) return;
+  finput = fopen(P->Call.ForcesFile, "r");
+        //fprintf(stderr, "File pointer says %p\n", finput);
+        if (finput == NULL)
+                Error(InitReading, "ForcesFile does not exist.");
+        fscanf(finput, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n", line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line);        // skip header line
+  for (i=0; i < I->Max_Types; i++) 
+    for (j=0; j < I->I[i].Max_IonsOfType; j++)
+                        if (feof(finput)) {
+                                Error(InitReading, "ForcesFile does not contain enough lines.");
+                        } else {
+                                //fscanf(finput, "%s\n", line);
+                                //if (strlen(line) < 100) {
+                                        //fprintf(stderr, "i %d, j %d, length of line %ld, line %s.\n", i,j, strlen(line), line);
+                                        //Error(InitReading, "ForcesFile does not contain enough lines or line is faulty.");
+                                //} else
+                                        //fprintf(stderr, "Parsed line: '%s'\n", line);
+                fscanf(finput, "%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n", &i1,&j1, 
+                      &R[0], &R[1], &R[2],
+                      &I->I[i].FIon[0+j*NDIM], &I->I[i].FIon[1+j*NDIM], &I->I[i].FIon[2+j*NDIM],
+                    &I->I[i].FIonL[0+j*NDIM], &I->I[i].FIonL[1+j*NDIM], &I->I[i].FIonL[2+j*NDIM],
+                  &I->I[i].FIonNL[0+j*NDIM], &I->I[i].FIonNL[1+j*NDIM], &I->I[i].FIonNL[2+j*NDIM],
+                &I->I[i].FMagnetic[0+j*NDIM], &I->I[i].FMagnetic[1+j*NDIM], &I->I[i].FMagnetic[2+j*NDIM],
+                      &I->I[i].FEwald[0+j*NDIM], &I->I[i].FEwald[1+j*NDIM], &I->I[i].FEwald[2+j*NDIM]);
+                                if ((i != (int)i1) || (j != (int)j1))
+                                        Error(InitReading, "Line does not match to desired ion.");
+                }
+        fclose(finput);
+  //fprintf(stderr, "MeanForce:\t%e\n", Run->MeanForce[0]);
+};
 
 /** Frees memory IonType.

pcp/src/ions.h

@@ -121 +121 @@
 /* Functions */
 void IonsInitRead(struct Problem *P, FILE *source);
+void InitThermostats(struct Problem *P, FILE *source);
 void CalculateEwald(struct Problem *P, int first);
 void RemoveIonsRead(struct Ions *I);