git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@332 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/min_cg.cpp

@@ -72,12 +72,12 @@ void MinCG::init()
   delete [] fixarg;
   fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
 
-  // zero local vectors before first atom exchange
+  // zero gradient vectors before first atom exchange
 
-  set_local_vectors();
+  setup_vectors();
   for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0;
 
-  // virial_thermo = how virial computed on thermo timesteps
+  // virial_thermo is how virial should be computed on thermo timesteps
   // 1 = computed explicity by pair, 2 = computed implicitly by pair
 
   if (force->newton_pair) virial_thermo = 2;
@@ -127,31 +127,28 @@ void MinCG::init()
 
 void MinCG::run()
 {
-  int i;
+  double tmp,*f;
-  double tmp;
 
   // set initial force & energy
 
   setup();
+  setup_vectors();
   output->thermo->compute_pe();
-  energy = output->thermo->potential_energy;
+  ecurrent = output->thermo->potential_energy;
-  energy += energy_extra;
 
   // stats for Finish to print
 	
-  einitial = energy;
+  einitial = ecurrent;
 
+  f = atom->f[0];
   tmp = 0.0;
-  for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
+  for (int i = 0; i < ndof; i++) tmp += f[i]*f[i];
   MPI_Allreduce(&tmp,&gnorm2_init,1,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < ndof_extra; i++) gnorm2_init += fextra[i]*fextra[i];
   gnorm2_init = sqrt(gnorm2_init);
 
   tmp = 0.0;
-  for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
+  for (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
   MPI_Allreduce(&tmp,&gnorminf_init,1,MPI_DOUBLE,MPI_MAX,world);
-  for (i = 0; i < ndof_extra; i++)
-    gnorminf_init = MAX(gnorminf_init,fabs(fextra[i]));
 
   // minimizer iterations
 
@@ -172,19 +169,16 @@ void MinCG::run()
     output->next_dump_any = update->ntimestep;
     if (output->restart_every) output->next_restart = update->ntimestep;
     output->next_thermo = update->ntimestep;
-    eng_force();
+    int ntmp;
+    double *xtmp,*htmp,etmp;
+    eng_force(&ntmp,&xtmp,&htmp,&etmp);
     output->write(update->ntimestep);
   }
   timer->barrier_stop(TIME_LOOP);
 
   // delete fix at end of run, so its atom arrays won't persist
-  // delete extra arrays
 
   modify->delete_fix("MINIMIZE");
-  delete [] xextra;
-  delete [] fextra;
-  delete [] gextra;
-  delete [] hextra;
 
   // reset reneighboring criteria
 
@@ -194,20 +188,17 @@ void MinCG::run()
 
   // stats for Finish to print
 	
-  efinal = energy;
+  efinal = ecurrent;
 
   f = atom->f[0];
   tmp = 0.0;
-  for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
+  for (int i = 0; i < ndof; i++) tmp += f[i]*f[i];
   MPI_Allreduce(&tmp,&gnorm2_final,1,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < ndof_extra; i++) gnorm2_final += fextra[i]*fextra[i];
   gnorm2_final = sqrt(gnorm2_final);
 
   tmp = 0.0;
-  for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
+  for (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
   MPI_Allreduce(&tmp,&gnorminf_final,1,MPI_DOUBLE,MPI_MAX,world);
-  for (i = 0; i < ndof_extra; i++)
-    gnorminf_final = MAX(gnorminf_final,fabs(fextra[i]));
 }
 
 /* ----------------------------------------------------------------------
@@ -216,20 +207,6 @@ void MinCG::run()
 
 void MinCG::setup()
 {
-  // allocate extra arrays
-  // couldn't do in init(), b/c modify and fixes weren't yet init()
-  // set initial xextra values via fixes
-
-  ndof_extra = modify->min_dof();
-  xextra = new double[ndof_extra];
-  fextra = new double[ndof_extra];
-  gextra = new double[ndof_extra];
-  hextra = new double[ndof_extra];
-
-  modify->min_xinitial(xextra);
-
-  // perform usual setup
-
   if (comm->me == 0 && screen) fprintf(screen,"Setting up minimization ...\n");
 
   // setup domain, communication and neighboring
@@ -245,7 +222,7 @@ void MinCG::setup()
   comm->borders();
   neighbor->build();
   neighbor->ncalls = 0;
-  set_local_vectors();
+  setup_vectors();
 
   // compute all forces
 
@@ -270,7 +247,6 @@ void MinCG::setup()
   if (force->newton) comm->reverse_communicate();
 
   modify->setup();
-  energy_extra = modify->min_energy(xextra,fextra);
   output->setup(1);
 }
 
@@ -283,14 +259,14 @@ void MinCG::iterate(int n)
 {
   int i,gradsearch,fail;
   double alpha,beta,gg,dot[2],dotall[2];
+  double *f;
 
-  for (i = 0; i < ndof; i++) h[i] = g[i] = f[i];
+  f = atom->f[0];
-  for (i = 0; i < ndof_extra; i++) hextra[i] = gextra[i] = fextra[i];
+  for (int i = 0; i < ndof; i++) h[i] = g[i] = f[i];
 
   dot[0] = 0.0;
   for (i = 0; i < ndof; i++) dot[0] += f[i]*f[i];
   MPI_Allreduce(dot,&gg,1,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < ndof_extra; i++) gg += fextra[i]*fextra[i];
 
   neval = 0;
   gradsearch = 1;
@@ -301,8 +277,8 @@ void MinCG::iterate(int n)
 
     // line minimization along direction h from current atom->x
 
-    eprevious = energy;
+    eprevious = ecurrent;
-    fail = (this->*linemin)(neval);
+    fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval);
 
     // if max_eval exceeded, all done
     // if linemin failed or energy did not decrease sufficiently:
@@ -311,8 +287,8 @@ void MinCG::iterate(int n)
 
     if (neval >= update->max_eval) break;
 
-    if (fail || fabs(energy-eprevious) <= 
+    if (fail || fabs(ecurrent-eprevious) <= 
-    	update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS)) {
+    	update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) {
       if (gradsearch == 1) break;
       gradsearch = -1;
     }
@@ -323,17 +299,14 @@ void MinCG::iterate(int n)
     // force new search dir to be grad dir if need to restart CG
     // set gradsearch to 1 if will search in grad dir on next iteration
 
+    f = atom->f[0];
     dot[0] = dot[1] = 0.0;
     for (i = 0; i < ndof; i++) {
       dot[0] += f[i]*f[i];
       dot[1] += f[i]*g[i];
     }
     MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
-    for (i = 0; i < ndof_extra; i++) {
+
-      dotall[0] += fextra[i]*fextra[i];
-      dotall[1] += fextra[i]*gextra[i];
-    }
-    
     beta = MAX(0.0,(dotall[0] - dotall[1])/gg);
     gg = dotall[0];
     if (gg < EPS) break;
@@ -346,10 +319,6 @@ void MinCG::iterate(int n)
       g[i] = f[i];
       h[i] = g[i] + beta*h[i];
     }
-    for (i = 0; i < ndof_extra; i++) {
-      gextra[i] = fextra[i];
-      hextra[i] = gextra[i] + beta*hextra[i];
-    }
 
     // output for thermo, dump, restart files
 
@@ -365,11 +334,9 @@ void MinCG::iterate(int n)
    set ndof and vector pointers after atoms have migrated
 ------------------------------------------------------------------------- */
 
-void MinCG::set_local_vectors()
+void MinCG::setup_vectors()
 {
   ndof = 3 * atom->nlocal;
-  x = atom->x[0];
-  f = atom->f[0];
   if (ndof) g = fix_minimize->gradient[0];
   else g = NULL;
   if (ndof) h = fix_minimize->searchdir[0];
@@ -379,11 +346,11 @@ void MinCG::set_local_vectors()
 /* ----------------------------------------------------------------------
    evaluate potential energy and forces
    may migrate atoms
-   energy = new objective function energy = poteng of atoms + eng_extra
+   new energy stored in ecurrent and returned (in case caller not in class)
-   atom->f, fextra = negative gradient of objective function
+   negative gradient will be stored in atom->f
 ------------------------------------------------------------------------- */
 
-void MinCG::eng_force()
+void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
 {
   // check for reneighboring
   // always communicate since minimizer moved atoms
@@ -408,7 +375,7 @@ void MinCG::eng_force()
     timer->stamp(TIME_COMM);
     neighbor->build();
     timer->stamp(TIME_NEIGHBOR);
-    set_local_vectors();
+    setup_vectors();
   }
 
   // eflag is always set, since minimizer needs potential energy
@@ -442,17 +409,21 @@ void MinCG::eng_force()
     timer->stamp(TIME_COMM);
   }
 
-  // min_post_force = forces on atoms that affect minimization
+  // fixes that affect minimization
-  // min_energy = energy, forces on extra degrees of freedom
 
   if (modify->n_min_post_force) modify->min_post_force(vflag);
-  if (modify->n_min_energy) energy_extra = modify->min_energy(xextra,fextra);
 
   // compute potential energy of system via Thermo
 
   output->thermo->compute_pe();
-  energy = output->thermo->potential_energy;
+  ecurrent = output->thermo->potential_energy;
-  energy += energy_extra;
+
+  // return updated ptrs to caller since atoms may have migrated
+
+  *pndof = ndof;
+  *px = atom->x[0];
+  *ph = h;
+  *peng = ecurrent;
 }
 
 /* ----------------------------------------------------------------------
@@ -516,15 +487,25 @@ void MinCG::force_clear(int vflag)
 
 /* ----------------------------------------------------------------------
    line minimization methods
-   find minimum-energy starting at x along h direction
+   find minimum-energy starting at x along dir direction
-   update atom->x by alpha, call eng_force() for result
+   input: n = # of degrees of freedom on this proc
-   alpha = distance moved along h to set x to minimun-energy configuration
+          x = ptr to atom->x[0] as vector
-   return 0 if successful move, 1 if failed (no move)
+	  dir = search direction as vector
-   insure last call to eng_force() is consistent with return
+	  eng = current energy at initial x
-     if fail, eng_force() of original x
+	  min/max dist = min/max distance to move any atom coord
-     if succeed, eng_force() at x + alpha*h
+   output: return 0 if successful move, set alpha
-   eng_force() may migrate atoms due to neighbor list build
+           return 1 if failed, no move, no need to set alpha
-     therefore linemin routines CANNOT store atom-based quantities
+           alpha = distance moved along dir to set x to min-eng config
+           caller has several quantities set via last call to eng_force()
+	     INSURE last call to eng_force() is consistent with returns
+	       if fail, eng_force() of original x
+	       if succeed, eng_force() at x + alpha*dir
+             atom->x = coords at new configuration
+	     atom->f = force (-Grad) is evaulated at new configuration
+	     ecurrent = energy of new configuration
+   NOTE: when call eng_force: n,x,dir,eng may change due to atom migration
+	 updated values are returned by eng_force()
+	 this routine CANNOT store atom-based quantities b/c of migration
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
@@ -534,55 +515,52 @@ void MinCG::force_clear(int vflag)
    quit as soon as energy starts to rise
 ------------------------------------------------------------------------- */
 
-int MinCG::linemin_scan(int &nfunc)
+int MinCG::linemin_scan(int n, double *x, double *dir, double eng,
+			double mindist, double maxdist,
+			double &alpha, int &nfunc)
 {
   int i;
-  double fmax,fme,elowest,alpha,alphamin,alphamax,alphalast;
+  double fmax,fme,elowest,alphamin,alphamax,alphalast;
 
   // alphamin = step that moves some atom coord by mindist
   // alphamax = step that moves some atom coord by maxdist
 
   fme = 0.0;
-  for (i = 0; i < ndof; i++) fme = MAX(fme,fabs(h[i]));
+  for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i]));
   MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world);
-  for (i = 0; i < ndof_extra; i++) fmax = MAX(fmax,fabs(hextra[i]));
-
   if (fmax == 0.0) return 1;
-  alphamin = dmin/fmax;
+
-  alphamax = dmax/fmax;
+  alphamin = mindist/fmax;
+  alphamax = maxdist/fmax;
 
   // if minstep is already uphill, fail
   // if eng increases, stop and return previous alpha
   // if alphamax, stop and return alphamax
 
-  elowest = energy;
+  elowest = eng;
   alpha = alphamin;
 
   while (1) {
-    for (i = 0; i < ndof; i++) x[i] += alpha*h[i];
+    for (i = 0; i < n; i++) x[i] += alpha*dir[i];
-    for (i = 0; i < ndof_extra; i++) xextra[i] += alpha*hextra[i];
+    eng_force(&n,&x,&dir,&eng);
-    eng_force();
     nfunc++;
 
-    if (alpha == alphamin && energy >= elowest) {
+    if (alpha == alphamin && eng >= elowest) {
-      for (i = 0; i < ndof; i++) x[i] -= alpha*h[i];
+      for (i = 0; i < n; i++) x[i] -= alpha*dir[i];
-      for (i = 0; i < ndof_extra; i++) xextra[i] -= alpha*hextra[i];
+      eng_force(&n,&x,&dir,&eng);
-      eng_force();
       nfunc++;
       return 1;
     }
-    if (energy > elowest) {
+    if (eng > elowest) {
-      for (i = 0; i < ndof; i++) x[i] += (alphalast-alpha)*h[i];
+      for (i = 0; i < n; i++) x[i] += (alphalast-alpha)*dir[i];
-      for (i = 0; i < ndof_extra; i++)
+      eng_force(&n,&x,&dir,&eng);
-	xextra[i] += (alphalast-alpha)*hextra[i];
-      eng_force();
       nfunc++;
       alpha = alphalast;
       return 0;
     }      
     if (alpha == alphamax) return 0;
 
-    elowest = energy;
+    elowest = eng;
     alphalast = alpha;
     alpha *= SCAN_FACTOR;
     if (alpha > alphamax) alpha = alphamax;
@@ -596,44 +574,43 @@ int MinCG::linemin_scan(int &nfunc)
    prevents successvive func evals further apart in x than maxdist
 ------------------------------------------------------------------------- */
 
-int MinCG::linemin_secant(int &nfunc)
+int MinCG::linemin_secant(int n, double *x, double *dir, double eng,
+			  double mindist, double maxdist,
+			  double &alpha, int &nfunc)
 {
   int i,iter;
-  double eta,eta_prev,sigma0,sigmamax,alpha,alphadelta,fme,fmax,dsq,e0,tmp;
+  double eta,eta_prev,sigma0,sigmamax,alphadelta,fme,fmax,dsq,e0,tmp;
+  double *f;
   double epssq = SECANT_EPS * SECANT_EPS;
 
   // stopping criterion for secant iterations
 
   fme = 0.0;
-  for (i = 0; i < ndof; i++) fme += h[i]*h[i];
+  for (i = 0; i < n; i++) fme += dir[i]*dir[i];
   MPI_Allreduce(&fme,&dsq,1,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < ndof_extra; i++) dsq += hextra[i]*hextra[i];
 
   // sigma0 = smallest allowed step of mindist
   // sigmamax = largest allowed step (in single iteration) of maxdist
 
   fme = 0.0;
-  for (i = 0; i < ndof; i++) fme = MAX(fme,fabs(h[i]));
+  for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i]));
   MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world);
-  for (i = 0; i < ndof_extra; i++) fmax = MAX(fmax,fabs(hextra[i]));
-
   if (fmax == 0.0) return 1;
-  sigma0 = dmin/fmax;
+
-  sigmamax = dmax/fmax;
+  sigma0 = mindist/fmax;
+  sigmamax = maxdist/fmax;
 
   // eval func at sigma0
   // test if minstep is already uphill
 
-  e0 = energy;
+  e0 = eng;
-  for (i = 0; i < ndof; i++) x[i] += sigma0*h[i];
+  for (i = 0; i < n; i++) x[i] += sigma0*dir[i];
-  for (i = 0; i < ndof_extra; i++) xextra[i] += sigma0*hextra[i];
+  eng_force(&n,&x,&dir,&eng);
-  eng_force();
   nfunc++;
 
-  if (energy >= e0) {
+  if (eng >= e0) {
-    for (i = 0; i < ndof; i++) x[i] -= sigma0*h[i];
+    for (i = 0; i < n; i++) x[i] -= sigma0*dir[i];
-    for (i = 0; i < ndof_extra; i++) xextra[i] -= sigma0*hextra[i];
+    eng_force(&n,&x,&dir,&eng);
-    eng_force();
     nfunc++;
     return 1;
   }
@@ -641,25 +618,24 @@ int MinCG::linemin_secant(int &nfunc)
   // secant iterations
   // alphadelta = new increment to move, alpha = accumulated move
 
+  f = atom->f[0];
   tmp = 0.0;
-  for (i = 0; i < ndof; i++) tmp -= f[i]*h[i];
+  for (i = 0; i < n; i++) tmp -= f[i]*dir[i];
   MPI_Allreduce(&tmp,&eta_prev,1,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < ndof_extra; i++) eta_prev -= fextra[i]*hextra[i];
 
   alpha = sigma0;
   alphadelta = -sigma0;
 
   for (iter = 0; iter < lineiter; iter++) {
     alpha += alphadelta;
-    for (i = 0; i < ndof; i++) x[i] += alphadelta*h[i];
+    for (i = 0; i < n; i++) x[i] += alphadelta*dir[i];
-    for (i = 0; i < ndof_extra; i++) xextra[i] += alphadelta*hextra[i];
+    eng_force(&n,&x,&dir,&eng);
-    eng_force();
     nfunc++;
 
+    f = atom->f[0];
     tmp = 0.0;
-    for (i = 0; i < ndof; i++) tmp -= f[i]*h[i];
+    for (i = 0; i < n; i++) tmp -= f[i]*dir[i];
     MPI_Allreduce(&tmp,&eta,1,MPI_DOUBLE,MPI_SUM,world);
-    for (i = 0; i < ndof_extra; i++) eta -= fextra[i]*hextra[i];
 
     alphadelta *= eta / (eta_prev - eta);
     eta_prev = eta;

src/min_cg.h

@@ -27,39 +27,32 @@ class MinCG : public Min {
   virtual void iterate(int);
 
  protected:
-  int virial_thermo;         // what vflag should be on thermo steps (1,2)
+  int virial_thermo;      // what vflag should be on thermo steps (1,2)
-  int pairflag,torqueflag,granflag;               // force clear flags
+  int pairflag,torqueflag,granflag;
   int neigh_every,neigh_delay,neigh_dist_check;   // copies of reneigh criteria
 
-  int maxpair;               // copies of Update quantities
+  int maxpair;            // copies of Update quantities
   double **f_pair;
 
   class FixMinimize *fix_minimize;  // fix that stores gradient vecs
-  double mindist,maxdist;    // min/max dist for coord delta in line search
+  double ecurrent;            // current potential energy
+  double mindist,maxdist;     // min/max dist for coord delta in line search
 
-  int ndof;                  // 3N degrees-of-freedom on this proc
+  int ndof;               // # of degrees-of-freedom on this proc
-  double *x;                 // vec of 3N coords, ptr to atom->x[0]
+  double *g,*h;           // local portion of gradient, searchdir vectors
-  double *f;                 // vec of 3N forces, ptr to atom->f[0]
-  double *g;                 // vec of 3N old forces, ptr to fix_minimize::g
-  double *h;                 // vec of 3N search dir, ptr to fix_minimize::h
 
-  int ndof_extra;            // extra degrees of freedom to include in min
+  typedef int (MinCG::*FnPtr)(int, double *, double *, double,
-  double energy_extra;       // extra potential energy
+			      double, double, double &, int &);
-  double *xextra;            // extra vecs associated with ndof_extra
+  FnPtr linemin;          // ptr to linemin functions
-  double *fextra;
-  double *gextra;
-  double *hextra;
 
-  double energy;             // potential energy of atoms and extra dof
+  int linemin_scan(int, double *, double *, double,
-
+		   double, double, double &, int &);
-  typedef int (MinCG::*FnPtr)(int &);
+  int linemin_secant(int, double *, double *, double,
-  FnPtr linemin;                        // ptr to linemin functions
+		     double, double, double &, int &);
-  int linemin_scan(int &);
-  int linemin_secant(int &);
 
   void setup();
-  void set_local_vectors();
+  void setup_vectors();
-  void eng_force();
+  void eng_force(int *, double **, double **, double *);
   void force_clear(int);
 };
 

src/min_cg_fr.cpp

@@ -39,14 +39,14 @@ void MinCGFR::iterate(int n)
 {
   int i,gradsearch,fail;
   double alpha,beta,gg,dot,dotall;
+  double *f;
 
+  f = atom->f[0];
   for (int i = 0; i < ndof; i++) h[i] = g[i] = f[i];
-  for (i = 0; i < ndof_extra; i++) hextra[i] = gextra[i] = fextra[i];
 
   dot = 0.0;
   for (i = 0; i < ndof; i++) dot += f[i]*f[i];
   MPI_Allreduce(&dot,&gg,1,MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < ndof_extra; i++) gg += fextra[i]*fextra[i];
 
   neval = 0;
   gradsearch = 1;
@@ -57,8 +57,8 @@ void MinCGFR::iterate(int n)
 
     // line minimization along direction h from current atom->x
 
-    eprevious = energy;
+    eprevious = ecurrent;
-    fail = (this->*linemin)(neval);
+    fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval);
 
     // if max_eval exceeded, all done
     // if linemin failed or energy did not decrease sufficiently:
@@ -67,8 +67,8 @@ void MinCGFR::iterate(int n)
 
     if (neval >= update->max_eval) break;
 
-    if (fail || fabs(energy-eprevious) <= 
+    if (fail || fabs(ecurrent-eprevious) <= 
-    	update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS)) {
+    	update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) {
       if (gradsearch == 1) break;
       gradsearch = -1;
     }
@@ -79,10 +79,10 @@ void MinCGFR::iterate(int n)
     // force new search dir to be grad dir if need to restart CG
     // set gradsesarch to 1 if will search in grad dir on next iteration
 
+    f = atom->f[0];
     dot = 0.0;
     for (i = 0; i < ndof; i++) dot += f[i]*f[i];
     MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
-    for (i = 0; i < ndof_extra; i++) dotall += fextra[i]*fextra[i];
 
     beta = dotall/gg;
     gg = dotall;
@@ -96,10 +96,6 @@ void MinCGFR::iterate(int n)
       g[i] = f[i];
       h[i] = g[i] + beta*h[i];
     }
-    for (i = 0; i < ndof_extra; i++) {
-      gextra[i] = fextra[i];
-      hextra[i] = gextra[i] + beta*hextra[i];
-    }
 
     // output for thermo, dump, restart files
 

src/min_sd.cpp

@@ -35,9 +35,10 @@ void MinSD::iterate(int n)
 {
   int i,fail;
   double alpha,dot,dotall;
+  double *f;
 
+  f = atom->f[0];
   for (int i = 0; i < ndof; i++) h[i] = f[i];
-  for (i = 0; i < ndof_extra; i++) hextra[i] = fextra[i];
 
   neval = 0;
 
@@ -47,29 +48,28 @@ void MinSD::iterate(int n)
 
     // line minimization along direction h from current atom->x
 
-    eprevious = energy;
+    eprevious = ecurrent;
-    fail = (this->*linemin)(neval);
+    fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval);
 
     // if max_eval exceeded, all done
     // if linemin failed or energy did not decrease sufficiently, all done
 
     if (neval >= update->max_eval) break;
 
-    if (fail || fabs(energy-eprevious) <= 
+    if (fail || fabs(ecurrent-eprevious) <= 
-    	update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS))
+    	update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS))
       break;
 
     // set h to new f = -Grad(x)
     // done if size sq of grad vector < EPS
 
+    f = atom->f[0];
     dot = 0.0;
     for (i = 0; i < ndof; i++) dot += f[i]*f[i];
     MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
-    for (i = 0; i < ndof_extra; i++) dotall += fextra[i]*fextra[i];
     if (dotall < EPS) break;
 
     for (i = 0; i < ndof; i++) h[i] = f[i];
-    for (i = 0; i < ndof_extra; i++) hextra[i] = fextra[i];
 
     // output for thermo, dump, restart files
 

src/style_class2.h

@@ -0,0 +1,56 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef AngleInclude
+#include "angle_class2.h"
+#endif
+
+#ifdef AngleClass
+AngleStyle(class2,AngleClass2)
+#endif
+
+#ifdef BondInclude
+#include "bond_class2.h"
+#endif
+
+#ifdef BondClass
+BondStyle(class2,BondClass2)
+#endif
+
+#ifdef DihedralInclude
+#include "dihedral_class2.h"
+#endif
+
+#ifdef DihedralClass
+DihedralStyle(class2,DihedralClass2)
+#endif
+
+#ifdef ImproperInclude
+#include "improper_class2.h"
+#endif
+
+#ifdef ImproperClass
+ImproperStyle(class2,ImproperClass2)
+#endif
+
+#ifdef PairInclude
+#include "pair_lj_class2.h"
+#include "pair_lj_class2_coul_cut.h"
+#include "pair_lj_class2_coul_long.h"
+#endif
+
+#ifdef PairClass
+PairStyle(lj/class2,PairLJClass2)
+PairStyle(lj/class2/coul/cut,PairLJClass2CoulCut)
+PairStyle(lj/class2/coul/long,PairLJClass2CoulLong)
+#endif

src/style_dpd.h

@@ -0,0 +1,28 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef AtomInclude
+#include "atom_vec_dpd.h"
+#endif
+
+#ifdef AtomClass
+AtomStyle(dpd,AtomVecDPD)
+#endif
+
+#ifdef PairInclude
+#include "pair_dpd.h"
+#endif
+
+#ifdef PairClass
+PairStyle(dpd,PairDPD)
+#endif

src/style_granular.h

@@ -0,0 +1,50 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef AtomInclude
+#include "atom_vec_granular.h"
+#endif
+
+#ifdef AtomClass
+AtomStyle(granular,AtomVecGranular)
+# endif
+
+#ifdef FixInclude
+#include "fix_freeze.h"
+#include "fix_gran_diag.h"
+#include "fix_nve_gran.h"
+#include "fix_pour.h"
+#include "fix_shear_history.h"
+#include "fix_wall_gran.h"
+#endif
+
+#ifdef FixClass
+FixStyle(freeze,FixFreeze)
+FixStyle(gran/diag,FixGranDiag)
+FixStyle(nve/gran,FixNVEGran)
+FixStyle(pour,FixPour)
+FixStyle(SHEAR_HISTORY,FixShearHistory)
+FixStyle(wall/gran,FixWallGran)
+#endif
+
+#ifdef PairInclude
+#include "pair_gran_hertzian.h"
+#include "pair_gran_history.h"
+#include "pair_gran_no_history.h"
+#endif
+
+#ifdef PairClass
+PairStyle(gran/hertzian,PairGranHertzian)
+PairStyle(gran/history,PairGranHistory)
+PairStyle(gran/no_history,PairGranNoHistory)
+#endif

src/style_kspace.h

@@ -0,0 +1,38 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef KSpaceInclude
+#include "ewald.h"
+#include "pppm.h"
+#include "pppm_tip4p.h"
+#endif
+
+#ifdef KSpaceClass
+KSpaceStyle(ewald,Ewald)
+KSpaceStyle(pppm,PPPM)
+KSpaceStyle(pppm/tip4p,PPPMTIP4P)
+#endif
+
+#ifdef PairInclude
+#include "pair_buck_coul_long.h"
+#include "pair_lj_cut_coul_long.h"
+#include "pair_lj_cut_coul_long_tip4p.h"
+#include "pair_lj_charmm_coul_long.h"
+#endif
+
+#ifdef PairClass
+PairStyle(buck/coul/long,PairBuckCoulLong)
+PairStyle(lj/cut/coul/long,PairLJCutCoulLong)
+PairStyle(lj/cut/coul/long/tip4p,PairLJCutCoulLongTIP4P)
+PairStyle(lj/charmm/coul/long,PairLJCharmmCoulLong)
+#endif

src/style_manybody.h

@@ -0,0 +1,28 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PairInclude
+#include "pair_eam.h"
+#include "pair_eam_alloy.h"
+#include "pair_eam_fs.h"
+#include "pair_sw.h"
+#include "pair_tersoff.h"
+#endif
+
+#ifdef PairClass
+PairStyle(eam,PairEAM)
+PairStyle(eam/alloy,PairEAMAlloy)
+PairStyle(eam/fs,PairEAMFS)
+PairStyle(sw,PairSW)
+PairStyle(tersoff,PairTersoff)
+#endif

src/style_molecule.h

@@ -0,0 +1,116 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef AngleInclude
+#include "angle_charmm.h"
+#include "angle_cosine.h"
+#include "angle_cosine_squared.h"
+#include "angle_harmonic.h"
+#include "angle_hybrid.h"
+#endif
+
+#ifdef AngleClass
+AngleStyle(charmm,AngleCharmm)
+AngleStyle(cosine,AngleCosine)
+AngleStyle(cosine/squared,AngleCosineSquared)
+AngleStyle(harmonic,AngleHarmonic)
+AngleStyle(hybrid,AngleHybrid)
+#endif
+
+#ifdef AtomInclude
+#include "atom_vec_angle.h"
+#include "atom_vec_bond.h"
+#include "atom_vec_full.h"
+#include "atom_vec_molecular.h"
+#endif
+
+#ifdef AtomClass
+AtomStyle(angle,AtomVecAngle)
+AtomStyle(bond,AtomVecBond)
+AtomStyle(full,AtomVecFull)
+AtomStyle(molecular,AtomVecMolecular)
+#endif
+
+#ifdef BondInclude
+#include "bond_fene.h"
+#include "bond_fene_expand.h"
+#include "bond_harmonic.h"
+#include "bond_hybrid.h"
+#include "bond_morse.h"
+#include "bond_nonlinear.h"
+#include "bond_quartic.h"
+#endif
+
+#ifdef BondClass
+BondStyle(fene,BondFENE)
+BondStyle(fene/expand,BondFENEExpand)
+BondStyle(harmonic,BondHarmonic)
+BondStyle(hybrid,BondHybrid)
+BondStyle(morse,BondMorse)
+BondStyle(nonlinear,BondNonlinear)
+BondStyle(quartic,BondQuartic)
+#endif
+
+#ifdef DihedralInclude
+#include "dihedral_charmm.h"
+#include "dihedral_harmonic.h"
+#include "dihedral_helix.h"
+#include "dihedral_hybrid.h"
+#include "dihedral_multi_harmonic.h"
+#include "dihedral_opls.h"
+#endif
+
+#ifdef DihedralClass
+DihedralStyle(charmm,DihedralCharmm)
+DihedralStyle(harmonic,DihedralHarmonic)
+DihedralStyle(helix,DihedralHelix)
+DihedralStyle(hybrid,DihedralHybrid)
+DihedralStyle(multi/harmonic,DihedralMultiHarmonic)
+DihedralStyle(opls,DihedralOPLS)
+#endif
+
+#ifdef DumpInclude
+#include "dump_bond.h"
+#endif
+
+#ifdef DumpClass
+DumpStyle(bond,DumpBond)
+#endif
+
+#ifdef FixInclude
+#endif
+
+#ifdef FixClass
+#endif
+
+#ifdef ImproperInclude
+#include "improper_cvff.h"
+#include "improper_harmonic.h"
+#include "improper_hybrid.h"
+#endif
+
+#ifdef ImproperClass
+ImproperStyle(cvff,ImproperCvff)
+ImproperStyle(harmonic,ImproperHarmonic)
+ImproperStyle(hybrid,ImproperHybrid)
+#endif
+
+#ifdef PairInclude
+#include "pair_lj_charmm_coul_charmm.h"
+#include "pair_lj_charmm_coul_charmm_implicit.h"
+#endif
+
+#ifdef PairClass
+PairStyle(lj/charmm/coul/charmm,PairLJCharmmCoulCharmm)
+PairStyle(lj/charmm/coul/charmm/implicit,PairLJCharmmCoulCharmmImplicit)
+#endif

src/style_opt.h

@@ -0,0 +1,30 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PairInclude
+#include "pair_eam_opt.h"
+#include "pair_eam_alloy_opt.h"
+#include "pair_eam_fs_opt.h"
+#include "pair_lj_charmm_coul_long_opt.h"
+#include "pair_lj_cut_opt.h"
+#include "pair_morse_opt.h"
+#endif
+
+#ifdef PairClass
+PairStyle(eam/opt,PairEAMOpt)
+PairStyle(eam/alloy/opt,PairEAMAlloyOpt)
+PairStyle(eam/fs/opt,PairEAMFSOpt)
+PairStyle(lj/cut/opt,PairLJCutOpt)
+PairStyle(lj/charmm/coul/long/opt,PairLJCharmmCoulLongOpt)
+PairStyle(morse/opt,PairMorseOpt)
+#endif

src/style_xtc.h

@@ -0,0 +1,20 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef DumpInclude
+#include "dump_xtc.h"
+#endif
+
+#ifdef DumpClass
+DumpStyle(xtc,DumpXTC)
+#endif