Changeset 5712cb for pcp/src


Ignore:
Timestamp:
Apr 21, 2008, 2:19:23 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
963310a
Parents:
5a538b
git-author:
Frederik Heber <heber@…> (04/18/08 14:01:58)
git-committer:
Frederik Heber <heber@…> (04/21/08 14:19:23)
Message:

CalculateForce(): added parsing of ForcesFile, dropped debug()s and PsiStep=MaxPsiStep, always InitDensity not doing StructOpt (perturbed change densities), commented-put TestSawtooth(), DoOutCurr only in uppermost level, print intermediate energy of not StructOpt/MD, do Force calculation at the very end

File:
1 edited

Legend:

Unmodified
Added
Removed

  • pcp/src/run.c

    r5a538b r5712cb  
    10731073  int Stop=0, SuperStop = 0, OuterStop = 0;
    10741074  //int i, j;
    1075   while ((R->LevSNo > R->InitLevSNo) || (!Stop && R->LevSNo == R->InitLevSNo)) {
    1076     // occupied
    1077     R->PsiStep = R->MaxPsiStep; // reset in-Psi-minimisation-counter, so that we really advance to the next wave function
    1078     R->OldActualLocalPsiNo = R->ActualLocalPsiNo; // reset OldActualLocalPsiNo, as it might still point to a perturbed wave function from last level
    1079     UpdateGramSchOldActualPsiNo(P,Psi);
    1080     MinimiseOccupied(P, &Stop, &SuperStop);
    1081     if (!I->StructOpt) {
    1082       if ((P->Call.ReadSrcFiles != 1) || (!ParseWannierFile(P))) { // only localize and store if they have just been minimised (hence don't come solely from file), otherwise read stored values from file (if possible)
    1083         SpeedMeasure(P, WannierTime, StartTimeDo);
    1084         ComputeMLWF(P);   // localization of orbitals
    1085         SpeedMeasure(P, WannierTime, StopTimeDo);
    1086         OutputVisSrcFiles(P, Occupied); // rewrite now localized orbitals
    1087   //      if (!TestReadnWriteSrcDensity(P,Occupied))
    1088   //        Error(SomeError,"TestReadnWriteSrcDensity failed!");
    1089       }
    1090    
     1075  SpeedMeasure(P, SimTime, StartTimeDo);
     1076  if ((F->DoOutVis == 2) || (P->Call.ForcesFile == NULL)) {  // if we want to draw those pretty density pictures, we have to solve the ground state in any case
     1077    while ((R->LevSNo > R->InitLevSNo) || (!Stop && R->LevSNo == R->InitLevSNo)) {
     1078      // occupied
     1079      //R->PsiStep = R->MaxPsiStep; // reset in-Psi-minimisation-counter, so that we really advance to the next wave function
     1080      R->OldActualLocalPsiNo = R->ActualLocalPsiNo; // reset OldActualLocalPsiNo, as it might still point to a perturbed wave function from last level
     1081      UpdateGramSchOldActualPsiNo(P,Psi);
     1082      ControlNativeDensity(P);
     1083      MinimiseOccupied(P, &Stop, &SuperStop);
     1084      if (!I->StructOpt) {
     1085        if ((P->Call.ReadSrcFiles != 1) || (!ParseWannierFile(P))) { // only localize and store if they have just been minimised (hence don't come solely from file), otherwise read stored values from file (if possible)
     1086          SpeedMeasure(P, WannierTime, StartTimeDo);
     1087          ComputeMLWF(P);   // localization of orbitals
     1088          SpeedMeasure(P, WannierTime, StopTimeDo);
     1089          OutputVisSrcFiles(P, Occupied); // rewrite now localized orbitals
     1090    //      if (!TestReadnWriteSrcDensity(P,Occupied))
     1091    //        Error(SomeError,"TestReadnWriteSrcDensity failed!");
     1092        }
     1093     
    10911094//      // plot psi cuts
    10921095//      for (i=0; i < Psi->MaxPsiOfType; i++)  // go through all wave functions (here without the extra ones for each process)
    … …  
    10981101//            }
    10991102         
    1100       // unoccupied calc
    1101       if (R->DoUnOccupied) {
    1102         MinimiseUnoccupied(P, &Stop, &SuperStop);
    1103       }
    1104       // hamiltonian
    1105       CalculateHamiltonian(P);  // lambda_{kl} needed (and for bandgap after UnOccupied)
    1106      
    1107       // perturbed calc
    1108       if ((R->DoPerturbation)) { // && R->LevSNo <= R->InitLevSNo) {
    1109         AllocCurrentDensity(R->Lev0->Dens);// lock current density arrays
    1110         MinimisePerturbed(P, &Stop, &SuperStop); // herein InitDensityCalculation() is called, thus no need to call it beforehand
    1111 
    1112         SpeedMeasure(P, CurrDensTime, StartTimeDo);
    1113         if (SuperStop != 1) {
    1114           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
    1115             R->DoFullCurrent = 1; // set to 1 if it was 2 but Check...() yielded necessity
    1116             debug(P,"Filling with Delta j ...");
    1117             FillDeltaCurrentDensity(P);
    1118           }// else
    1119             //debug(P,"There is no overlap between orbitals.");
    1120           //debug(P,"Filling with j ...");
    1121           //FillCurrentDensity(P);
     1103        // unoccupied calc
     1104        if (R->DoUnOccupied) {
     1105          MinimiseUnoccupied(P, &Stop, &SuperStop);
    11221106        }
    1123         SpeedMeasure(P, CurrDensTime, StopTimeDo);
    1124         TestCurrent(P,0);
    1125         TestCurrent(P,1);
    1126         TestCurrent(P,2);
    1127         if (F->DoOutCurr) {
    1128           debug(P,"OutputCurrentDensity");
    1129           OutputCurrentDensity(P);
    1130         }
    1131         if (R->VectorPlane != -1) {
    1132           debug(P,"PlotVectorPlane");
    1133           PlotVectorPlane(P,R->VectorPlane,R->VectorCut);
    1134         }
    1135         fprintf(stderr,"(%i) ECut [L%i] = %e Ht\n", P->Par.me, R->Lev0->LevelNo, pow(2*M_PI*M_PI/Lat->Volume*R->Lev0->TotalAllMaxG, 2./3.));
    1136         debug(P,"Calculation of magnetic susceptibility");
    1137         CalculateMagneticSusceptibility(P);
    1138         debug(P,"Normal calculation of shielding over R-space");
    1139         CalculateChemicalShielding(P);
    1140         debug(P,"Reciprocal calculation of shielding over G-space");
    1141         CalculateChemicalShieldingByReciprocalCurrentDensity(P);
    1142         SpeedMeasure(P, CurrDensTime, StopTimeDo);
    1143         DisAllocCurrentDensity(R->Lev0->Dens);    // unlock current density arrays
    1144       } else {
     1107        // hamiltonian
     1108        CalculateHamiltonian(P);  // lambda_{kl} needed (and for bandgap after UnOccupied)
     1109
     1110        //TestSawtooth(P, 0);
     1111        //TestSawtooth(P, 1);
     1112        //TestSawtooth(P, 2);
     1113               
     1114        // perturbed calc
     1115        if ((R->DoPerturbation)) { // && R->LevSNo <= R->InitLevSNo) {
     1116          AllocCurrentDensity(R->Lev0->Dens);// lock current density arrays
     1117          MinimisePerturbed(P, &Stop, &SuperStop); // herein InitDensityCalculation() is called, thus no need to call it beforehand
     1118
     1119          SpeedMeasure(P, CurrDensTime, StartTimeDo);
     1120          if (SuperStop != 1) {
     1121            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
     1122              R->DoFullCurrent = 1; // set to 1 if it was 2 but Check...() yielded necessity
     1123              debug(P,"Filling with Delta j ...");
     1124              FillDeltaCurrentDensity(P);
     1125            }// else
     1126              //debug(P,"There is no overlap between orbitals.");
     1127            //debug(P,"Filling with j ...");
     1128            //FillCurrentDensity(P);
     1129          }
     1130          SpeedMeasure(P, CurrDensTime, StopTimeDo);
     1131          TestCurrent(P,0);
     1132          TestCurrent(P,1);
     1133          TestCurrent(P,2);
     1134          if (F->DoOutCurr) {
     1135            debug(P,"OutputCurrentDensity");
     1136            OutputCurrentDensity(P);
     1137          }
     1138          if (R->VectorPlane != -1) {
     1139            debug(P,"PlotVectorPlane");
     1140            PlotVectorPlane(P,R->VectorPlane,R->VectorCut);
     1141          }
     1142          fprintf(stderr,"(%i) ECut [L%i] = %e Ht\n", P->Par.me, R->Lev0->LevelNo, pow(2*M_PI*M_PI/Lat->Volume*R->Lev0->TotalAllMaxG, 2./3.));
     1143          CalculateMagneticSusceptibility(P);
     1144          debug(P,"Normal calculation of shielding over R-space");
     1145          CalculateChemicalShielding(P);
     1146          CalculateChemicalShieldingByReciprocalCurrentDensity(P);
     1147          SpeedMeasure(P, CurrDensTime, StopTimeDo);
     1148          DisAllocCurrentDensity(R->Lev0->Dens);    // unlock current density arrays
     1149        }  // end of if perturbation
    11451150        InitDensityCalculation(P);  // all unperturbed(!) wave functions've "changed" from ComputeMLWF(), thus reinit density
    1146       }
    1147       //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);
     1151      } else  // end of if StructOpt or MaxOuterStep
     1152                OutputVisSrcFiles(P, Occupied); // in structopt or MD write for every level
     1153
     1154      if ((!I->StructOpt) && (!R->MaxOuterStep))  // print intermediate levels energy results if we don't do MD or StructOpt
     1155              EnergyOutput(P, 1);
     1156      // next level
     1157      ChangeToLevUp(P, &Stop);
     1158      //if (isnan(LevS->LPsi->LocalPsi[R->ActualLocalPsiNo][0].re)) { fprintf(stderr,"(%i) WARNING in ChangeToLevUp(): LPsi->LocalPsi[%i]_[%i] = NaN!\n", P->Par.me, R->ActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }     
     1159      LevS = R->LevS; // re-set pointer that's changed from LevUp
    11481160    }
    1149      
    1150 //    if (!I->StructOpt && R->DoPerturbation) {
    1151 //      InitDensityCalculation(P);  // most of the density array were used during FillCurrentDensity(), thus reinit density
    1152 //    }
    1153 
    1154     // next level
    1155     ChangeToLevUp(P, &Stop);
    1156     //if (isnan(LevS->LPsi->LocalPsi[R->ActualLocalPsiNo][0].re)) { fprintf(stderr,"(%i) WARNING in ChangeToLevUp(): LPsi->LocalPsi[%i]_[%i] = NaN!\n", P->Par.me, R->ActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }     
    1157     LevS = R->LevS; // re-set pointer that's changed from LevUp
    1158   }
    1159   //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);
    1160   // necessary for correct ionic forces ...
    1161   SpeedMeasure(P, LocFTime, StartTimeDo);
    1162   CalculateIonLocalForce(P);
    1163   SpeedMeasure(P, LocFTime, StopTimeDo);
    1164   SpeedMeasure(P, NonLocFTime, StartTimeDo);
    1165   CalculateIonNonLocalForce(P);
    1166   SpeedMeasure(P, NonLocFTime, StopTimeDo);
    1167   CalculateEwald(P, 0);
    1168   CalculateIonForce(P);
     1161          SpeedMeasure(P, SimTime, StopTimeDo);
     1162          ControlNativeDensity(P);
     1163          TestGramSch(P,LevS,Psi, Occupied);   
     1164    // necessary for correct ionic forces ...
     1165    SpeedMeasure(P, LocFTime, StartTimeDo);
     1166    CalculateIonLocalForce(P);
     1167    SpeedMeasure(P, LocFTime, StopTimeDo);
     1168    SpeedMeasure(P, NonLocFTime, StartTimeDo);
     1169    CalculateIonNonLocalForce(P);
     1170    SpeedMeasure(P, NonLocFTime, StopTimeDo);
     1171    CalculateEwald(P, 1);
     1172    CalculateIonForce(P);
     1173  }
     1174  if (P->Call.ForcesFile != NULL) { // if we parse forces from file, values are written over (in case of DoOutVis)
     1175    fprintf(stderr, "Parsing Forces from file.\n");
     1176    //ParseIonForce(P);
     1177    //CalculateIonForce(P);
     1178  }
    11691179  CorrectForces(P);
    11701180  // ... on output of densities
    … …  
    11741184  }
    11751185
    1176   OutputNorm(stderr, P);
     1186  //OutputNorm(P);
    11771187  //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);
    1178   OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
    1179   TestGramSch(P,LevS,Psi, -1);
    1180   SpeedMeasure(P, SimTime, StopTimeDo);
    1181   /*TestGramSch(P, R->LevS, &P->Lat.Psi); */
    1182   ControlNativeDensity(P);
    1183   SpeedMeasure(P, LocFTime, StartTimeDo);
    1184   CalculateIonLocalForce(P);
    1185   SpeedMeasure(P, LocFTime, StopTimeDo);
    1186   SpeedMeasure(P, NonLocFTime, StartTimeDo);
    1187   CalculateIonNonLocalForce(P);
    1188   SpeedMeasure(P, NonLocFTime, StopTimeDo);
    1189   CalculateEwald(P, 0);
    1190   CalculateIonForce(P);
    1191   CorrectForces(P);
     1188  //OutputVis(P, P->R.Lev0->Dens->DensityArray[TotalDensity]);
     1189  /*TestGramSch(P,LevS,Psi, &P->Lat.Psi); */
    11921190  GetOuterStop(P);
    11931191  //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);
Note: See TracChangeset for help on using the changeset viewer.