Changeset c5bdb23 for pcp


Ignore:
Timestamp:
Apr 21, 2008, 2:19:25 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
69eca8
Parents:
3716b2
git-author:
Frederik Heber <heber@…> (04/18/08 16:14:42)
git-committer:
Frederik Heber <heber@…> (04/21/08 14:19:25)
Message:

CalculateEwald(): corrected force
RemoveConstraintItem() added

File:
1 edited

Legend:

Unmodified
Added
Removed

  • pcp/src/ions.c

    r3716b2 rc5bdb23  
    438438
    439439/** Calculated Ewald force.
     440 * The basic idea of the ewald method is to add a countercharge - here of gaussian form -
     441 * which shields the long-range charges making them effectively short-ranged, while the
     442 * countercharges themselves can be analytically (here by (hidden) Fouriertransform: it's
     443 * added to the electron density) integrated and withdrawn.
    440444 * \f[
    441445 *  E_{K-K} = \frac{1}{2} \sum_L \sum_{K=\{(I_s,I_a,J_s,J_a)|R_{I_s,I_a} - R_{J_s,J_a} - L\neq 0\}} \frac{Z_{I_s} Z_{J_s}} {|R_{I_s,I_a} - R_{J_s,J_a} - L|}
    … …  
    448452 * in a circular motion from the current cell to the outside up to Ions::R_cut.
    449453 * \param *P Problem at hand
    450  * \param first additional calculation beforehand if != 0
     454 * \param first if != 0 table mirror cells to take into (local!) account of ewald summation
    451455 */
    452456void CalculateEwald(struct Problem *P, int first)
    453457{
    454   int r,is,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec;
     458  int r,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec; //,is
    455459  int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
    456460  double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
    … …  
    462466  struct PseudoPot *PP = &P->PP;
    463467  struct Ions *I = &P->Ion;
    464   if (I->R_cut != 0.0) {
    465     if (first) {
    466       SpeedMeasure(P, InitSimTime, StopTimeDo);
    467       rcut2 = I->R_cut*I->R_cut;
    468       ActualVec =0;
    469       MaxCell = (int)5.*I->R_cut/pow(L->Volume,1./3.);
    470       if (MaxCell < 2) MaxCell = 2;
    471       for (i = -MaxCell; i <= MaxCell; i++) {
    472         for (j = -MaxCell; j <= MaxCell; j++) {
    473                                 for (k = -MaxCell; k <= MaxCell; k++) {
    474                                   r2 = 0.0;
    475                                   for (ir=0; ir <3; ir++) {
    476                                     t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
    477                                     r2 = t[ir]*t[ir];
    478                                   }
    479                                   Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
    480                                   if ((r2 <= rcut2) || Is_Neighb_Cell) {
    481                                     ActualVec++;
    482                                   }
    483                                 }
    484         }
    485       }
    486       MaxVec = ActualVec;
    487       I->MaxVec = ActualVec;
    488       MaxLocalVec = MaxVec / MaxPar;
    489       StartVec = ParMe*MaxLocalVec;
    490       r = MaxVec % MaxPar;
    491       if (ParMe >= r) {
    492         StartVec += r;
    493       } else {
    494         StartVec += ParMe;
    495         MaxLocalVec++;
    496       }
    497       I->MaxLocalVec = MaxLocalVec;
    498       LocalActualVec = ActualVec = 0;
     468  if (first) {
     469    SpeedMeasure(P, InitSimTime, StopTimeDo);
     470    rcut2 = I->R_cut*I->R_cut;
     471    ActualVec = 0;  // counts number of cells taken into account
     472    MaxCell = (int)(5.*I->R_cut/pow(L->Volume,1./3.));    // the number of cells in each direction is prop. to axis length over cutoff
     473    if (MaxCell < 2) MaxCell = 2;     // take at least 2
     474    for (i = -MaxCell; i <= MaxCell; i++) {
     475      for (j = -MaxCell; j <= MaxCell; j++) {
     476                                for (k = -MaxCell; k <= MaxCell; k++) {
     477                                  r2 = 0.;
     478                                  for (ir=0; ir <NDIM; ir++) {  // t is the offset vector pointing to distant cell (only diagonal entries)
     479                                    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
     480                                    r2 = t[ir]*t[ir];   // is squared distance to other cell
     481                                  }
     482                                  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));   // whether it's directly adjacent
     483                                  if ((r2 <= rcut2) || Is_Neighb_Cell) {    // if it's either adjacent or within cutoff, count it in
     484                                    ActualVec++;
     485                                  }
     486                                }
     487      }
     488    }
     489    MaxVec = ActualVec;   // store number of cells
     490    I->MaxVec = ActualVec;
     491    MaxLocalVec = MaxVec / MaxPar;  // share among processes
     492    StartVec = ParMe*MaxLocalVec; // for this process start at its patch
     493    r = MaxVec % MaxPar;  // rest not fitting...
     494    if (ParMe >= r) {   // ones who don't get extra cells have to shift their start vector
     495      StartVec += r;
     496    } else {    // others get extra cell and have to shift also by a bit
     497      StartVec += ParMe;
     498      MaxLocalVec++;
     499    }
     500    I->MaxLocalVec = MaxLocalVec;
     501    LocalActualVec = ActualVec = 0;
     502    // now go through the same loop again and note down every offset vector for cells in this process' patch
     503    if (MaxLocalVec != 0) {
    499504      I->RLatticeVec = (double *)
    500         Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald:");
    501       for (i = -MaxCell; i <= MaxCell; i++) {
    502         for (j = -MaxCell; j <= MaxCell; j++) {
    503                                 for (k = -MaxCell; k <= MaxCell; k++) {
    504                                   r2 = 0.0;
    505                                   for (ir=0; ir <3; ir++) {
    506                                     t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
    507                                     r2 = t[ir]*t[ir];
    508                                   }
    509                                   Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
    510                                   if ((r2 <= rcut2) || Is_Neighb_Cell) {
    511                                     if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) {
    512                                       I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
    513                                       I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
    514                                       I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
    515                                       LocalActualVec++;
    516                                     }
    517                                     ActualVec++;
    518                                   }
    519                                 }
    520         }
    521       }
    522       SpeedMeasure(P, InitSimTime, StartTimeDo);
    523     }
    524     SpeedMeasure(P, EwaldTime, StartTimeDo);
    525     esr = 0.0;
    526     for (is1=0;is1 < I->Max_Types; is1++)
    527       for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
    528         for (i=0; i< NDIM; i++)
    529                                 I->FTemp[i+NDIM*(ia1)] = 0;
    530     for (is1=0;is1 < I->Max_Types; is1++) {
    531       for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
    532         for (i=0;i<3;i++) {
    533                                 R1[i]=I->I[is1].R[i+NDIM*ia1];
    534         }
    535         for (ir2=0;ir2<I->MaxLocalVec;ir2++) {
    536                                 for (is2=0;is2 < I->Max_Types; is2++) {
    537                                   gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
    538                                   for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
    539                                     for (i=0;i<3;i++) {
    540                                       R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2];
    541                                     }
    542                                     erre2=0.0;
    543                                     for (i=0;i<3;i++) {
    544                                       R3[i] = R1[i]-R2[i];
    545                                       erre2+=R3[i]*R3[i];
    546                                     }
    547                                     if (erre2 > MYEPSILON) {
    548                                       arg=sqrt(erre2);
    549                                       fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;
    550                                      
    551                                       arg *= gkl;
    552                                       addesr = derf(arg);
    553                                       addesr = (1.0-addesr)*fac;
    554                                       esr += addesr;
    555                                       addpre=exp(-arg*arg)*gkl;
    556                                       addpre=PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;
    557                                       repand=(addesr+addpre)/erre2;
    558                                       for (i=0;i<3;i++) {
    559                                                                 fxx=repand*R3[i];
    560                                                                 /*fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);*/
    561                                                                 I->FTemp[i+NDIM*(ia1)] += fxx;
    562                                                                 I->FTemp[i+NDIM*(ia2)] -= fxx;
    563                                       }
    564                                     }
    565                                   }
    566                                 }
    567         }
    568       }
    569     }
    570   } else {
    571     esr = 0.0;
    572     for (is1=0;is1 < I->Max_Types; is1++)
    573       for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
    574         for (i=0; i< NDIM; i++)
    575           I->FTemp[i+NDIM*(ia1)] = 0.;
     505        Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald: I->RLatticeVec");
     506    } else {
     507      fprintf(stderr,"(%i) Warning in Ewald summation: Got MaxLocalVec == 0\n", P->Par.me);
     508    }
     509    for (i = -MaxCell; i <= MaxCell; i++) {
     510      for (j = -MaxCell; j <= MaxCell; j++) {
     511                                for (k = -MaxCell; k <= MaxCell; k++) {
     512                                  r2 = 0.;
     513                                  for (ir=0; ir <NDIM; ir++) {  // create offset vector from diagonal entries
     514                                    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
     515                                    r2 = t[ir]*t[ir];
     516                                  }
     517                                  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
     518                                  if ((r2 <= rcut2) || Is_Neighb_Cell) {  // if adjacent or within cutoff ...
     519                                    if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) { // ... and belongs to patch, store 'em
     520                                      I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
     521                                      I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
     522                                      I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
     523                                      LocalActualVec++;
     524                                    }
     525                                    ActualVec++;
     526                                  }
     527                                }
     528      }
     529    }
     530    SpeedMeasure(P, InitSimTime, StartTimeDo);
     531  }
     532  SpeedMeasure(P, EwaldTime, StartTimeDo);
     533
     534  // first set everything to zero
     535  esr = 0.0;
     536  //for (is1=0;is1 < I->Max_Types; is1++)
     537  // then take each ion of each type once
     538  for (is1=0;is1 < I->Max_Types; is1++) {
     539    // reset FTemp for IonType
     540    for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
     541      for (i=0; i< NDIM; i++)
     542        I->FTemp[i+NDIM*(ia1)] = 0.;
     543    // then sum for each on of the IonType
     544    for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
     545      for (i=0;i<NDIM;i++) {
     546                                R1[i]=I->I[is1].R[i+NDIM*ia1];    // R1 is local coordinate of current first ion
     547      }
     548      for (ir2=0;ir2<I->MaxLocalVec;ir2++) {  // for every cell of local patch (each with the same atoms)
     549                                for (is2=0;is2 < I->Max_Types; is2++) { // and for each ion of each type a second time
     550          // gkl = 1./(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
     551                                  gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
     552                                  for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
     553                                    for (i=0;i<3;i++) {
     554                                      R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2]; // R2 is coordinate of current second ion in the distant cell!
     555                                    }
     556
     557                                    erre2=0.0;
     558                                    for (i=0;i<NDIM;i++) {
     559                                      R3[i] = R1[i]-R2[i];    // R3 = R_is1,ia1 - R_is2,ia2 - L
     560                                      erre2+=R3[i]*R3[i];     // erre2 = |R_is1,ia1 - R_is2,ia2 - L|^2
     561                                    }
     562                                    if (erre2 > MYEPSILON) {  // appears as one over ... in fac, thus check if zero
     563              arg = sqrt(erre2); // arg = |R_is1,ia1 - R_is2,ia2 - L|
     564                                      fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;    // fac = -1/2 * Z_is1 * Z_is2 / arg
     565                                     
     566                                      arg *= gkl;             // arg = |R_is1,ia1 - R_is2,ia2 - L|/(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
     567                                      addesr = derf(arg); 
     568              addesr = (1.-addesr)*fac; // complementary error function: addesr = 1/2*erfc(arg)/|R_is1,ia1 - R_is2,ia2 - L|
     569                                      esr += addesr;
     570              addpre=exp(-arg*arg)*gkl;
     571              addpre = PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;  // addpre = exp(-|R_is1,ia1 - R_is2,ia2 - L|^2/(r_is1,gauss^2 + r_is2,gauss^2))/(sqrt(pi)*sqrt(r_is1,gauss^2 + r_is2,gauss^2))
     572              repand = (addesr+addpre)/erre2;
     573//              if ((fabs(repand) > MYEPSILON)) {
     574//                fprintf(stderr,"(%i) Ewald[is%d,ia%d/is%d,ia%d,(%d)]: Z1 %lg, Z2 %lg\tpre %lg\t esr %lg\t fac %lg\t arg %lg\tR3 (%lg,%lg,%lg)\n", P->Par.me,  is1, ia1, is2, ia2, ir2, PP->zval[is1], PP->zval[is2],addpre,addesr, fac, arg, R3[0], R3[1], R3[2]);
     575//              }
     576                                      for (i=0;i<NDIM;i++) {
     577                                                                fxx=repand*R3[i];
     578//                if ((fabs(repand) > MYEPSILON)) {
     579//                                                                fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);
     580//                }
     581                                                                I->FTemp[i+NDIM*(ia1)] += fxx*2.0;      // actio = reactio?
     582                                                                //I->FTemp[i+NDIM*(ia2)] -= fxx;
     583                                      }
     584                                    }
     585                                  }
     586                                }
     587      }
     588    }
     589    MPI_Allreduce (I->FTemp , I->I[is1].FEwald, NDIM*I->I[is1].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm); 
    576590  }
    577591  SpeedMeasure(P, EwaldTime, StopTimeDo);
    578   for (is=0;is < I->Max_Types; is++) {
    579     MPI_Allreduce (I->FTemp , I->I[is].FEwald, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
    580   }
     592  //for (is=0;is < I->Max_Types; is++) {
     593  //}
    581594  MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
     595//  fprintf(stderr,"\n");
    582596}
    583597
Note: See TracChangeset for help on using the changeset viewer.