From 1445a6d5368e8a3250432d0892e4a7429290bf46 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 21 Feb 2007 20:57:31 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@331 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_variable_atom.cpp | 96 ++++++++ src/compute_variable_atom.h | 37 +++ src/fix_ave_spatial.cpp | 439 ++++++++++++++++++++++++++++++++++ src/fix_ave_spatial.h | 52 ++++ src/fix_ave_time.cpp | 166 +++++++++++++ src/fix_ave_time.h | 46 ++++ src/min_cg.cpp | 222 +++++++++-------- src/min_cg.h | 39 +-- src/min_cg_fr.cpp | 18 +- src/min_sd.cpp | 14 +- src/style_class2.h | 56 ----- src/style_dpd.h | 28 --- src/style_granular.h | 50 ---- src/style_kspace.h | 38 --- src/style_manybody.h | 28 --- src/style_molecule.h | 116 --------- src/style_opt.h | 30 --- src/style_xtc.h | 20 -- 18 files changed, 1000 insertions(+), 495 deletions(-) create mode 100644 src/compute_variable_atom.cpp create mode 100644 src/compute_variable_atom.h create mode 100644 src/fix_ave_spatial.cpp create mode 100644 src/fix_ave_spatial.h create mode 100644 src/fix_ave_time.cpp create mode 100644 src/fix_ave_time.h diff --git a/src/compute_variable_atom.cpp b/src/compute_variable_atom.cpp new file mode 100644 index 0000000000..2db88251f4 --- /dev/null +++ b/src/compute_variable_atom.cpp @@ -0,0 +1,96 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_variable_atom.h" +#include "atom.h" +#include "input.h" +#include "variable.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeVariableAtom::ComputeVariableAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all("Illegal compute variable/atom command"); + + // store variable name + + int n = strlen(arg[3]) + 1; + varname = new char[n]; + strcpy(varname,arg[3]); + + peratom_flag = 1; + size_peratom = 0; + + nmax = 0; + result = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeVariableAtom::~ComputeVariableAtom() +{ + delete [] varname; + memory->sfree(result); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVariableAtom::init() +{ + // set ivariable used by this compute + + ivariable = input->variable->find(varname); + if (ivariable < 0) + error->all("Could not find compute variable name"); + + // test if variable of correct ATOM type +} + +/* ---------------------------------------------------------------------- */ + +void ComputeVariableAtom::compute_peratom() +{ + // grow result array if necessary + + if (atom->nlocal > nmax) { + memory->sfree(result); + nmax = atom->nmax; + result = (double *) memory->smalloc(nmax*sizeof(double), + "compute/variable/atom:result"); + scalar_atom = result; + } + + // parse variable once to create parse tree + // evaluate tree for all atoms, will be zero for atoms not in group + // free parse tree memory stored by Variable + + input->variable->build_parse_tree(ivariable); + input->variable->evaluate_parse_tree(igroup,result); + input->variable->free_parse_tree(); +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +int ComputeVariableAtom::memory_usage() +{ + int bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/compute_variable_atom.h b/src/compute_variable_atom.h new file mode 100644 index 0000000000..b290dc9846 --- /dev/null +++ b/src/compute_variable_atom.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_VARIABLE_ATOM_H +#define COMPUTE_VARIABLE_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeVariableAtom : public Compute { + public: + ComputeVariableAtom(class LAMMPS *, int, char **); + ~ComputeVariableAtom(); + void init(); + void compute_peratom(); + int memory_usage(); + + private: + int nmax,ivariable; + char *varname; + double *result; +}; + +} + +#endif diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp new file mode 100644 index 0000000000..69523e8789 --- /dev/null +++ b/src/fix_ave_spatial.cpp @@ -0,0 +1,439 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Pieter in't Veld (SNL) +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "fix_ave_spatial.h" +#include "atom.h" +#include "update.h" +#include "domain.h" +#include "lattice.h" +#include "modify.h" +#include "compute.h" +#include "group.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{LOWER,CENTER,UPPER,COORD}; +enum{DENSITY_MASS,DENSITY_NUM,VX,VY,VZ,FX,FY,FZ,COMPUTE}; +enum{SAMPLE,ALL}; + +/* ---------------------------------------------------------------------- */ + +FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 11) error->all("Illegal fix ave/spatial command"); + + nevery = atoi(arg[3]); + nfreq = atoi(arg[4]); + + if (strcmp(arg[5],"x") == 0) dim = 0; + else if (strcmp(arg[5],"y") == 0) dim = 1; + else if (strcmp(arg[5],"z") == 0) dim = 2; + else error->all("Illegal fix ave/spatial command"); + + if (strcmp(arg[6],"lower") == 0) originflag = LOWER; + if (strcmp(arg[6],"center") == 0) originflag = CENTER; + if (strcmp(arg[6],"upper") == 0) originflag = UPPER; + else originflag = COORD; + if (originflag == COORD) origin = atof(arg[6]); + + delta = atof(arg[7]); + + MPI_Comm_rank(world,&me); + if (me == 0) { + fp = fopen(arg[8],"w"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix ave/spatial file %s",arg[8]); + error->one(str); + } + } + + if (strcmp(arg[9],"density") == 0) { + if (strcmp(arg[10],"mass") == 0) which = DENSITY_MASS; + else if (strcmp(arg[10],"number") == 0) which = DENSITY_NUM; + else error->all("Illegal fix ave/spatial command"); + + } else if (strcmp(arg[9],"atom") == 0) { + if (strcmp(arg[10],"vx") == 0) which = VX; + else if (strcmp(arg[10],"vy") == 0) which = VY; + else if (strcmp(arg[10],"vz") == 0) which = VZ; + else if (strcmp(arg[10],"fx") == 0) which = FX; + else if (strcmp(arg[10],"fy") == 0) which = FY; + else if (strcmp(arg[10],"fz") == 0) which = FZ; + else error->all("Illegal fix ave/spatial command"); + + } else if (strcmp(arg[9],"compute") == 0) { + which = COMPUTE; + int n = strlen(arg[10]) + 1; + id_compute = new char[n]; + strcpy(id_compute,arg[10]); + + } else error->all("Illegal fix ave/spatial command"); + + // parse optional args + + normflag = ALL; + int scaleflag = 1; + + int iarg = 11; + while (iarg < narg) { + if (strcmp(arg[iarg],"norm") == 0) { + if (iarg+2 > narg) error->all("Illegal fix ave/spatial command"); + if (strcmp(arg[iarg+1],"all") == 0) normflag = ALL; + else if (strcmp(arg[iarg+1],"sample") == 0) normflag = SAMPLE; + else error->all("Illegal fix ave/spatial command"); + iarg += 2; + } else if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all("Illegal fix ave/spatial command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all("Illegal fix ave/spatial command"); + iarg += 2; + } else error->all("Illegal fix ave/spatial command"); + } + + // if density, no normalization by atom count should be done + // thus ALL and SAMPLE should give same answer, but code does normalize + // thus only ALL is computed correctly, so force norm to be ALL + + if (which == DENSITY_MASS || which == DENSITY_NUM) normflag = ALL; + + // setup scaling + + if (scaleflag && domain->lattice == NULL) + error->all("Use of fix ave/spatial with undefined lattice"); + + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + + // apply scaling factors + + double scale; + if (dim == 0) scale = xscale; + if (dim == 1) scale = yscale; + if (dim == 2) scale = zscale; + delta *= scale; + if (originflag == COORD) origin *= scale; + + // setup and error check + + if (nevery <= 0) error->all("Illegal fix ave/spatial command"); + if (nfreq < nevery || nfreq % nevery) + error->all("Illegal fix ave/spatial command"); + + if (delta <= 0.0) error->all("Illegal fix ave/spatial command"); + invdelta = 1.0/delta; + + nvalues = 1; + if (which == COMPUTE) { + int icompute = modify->find_compute(id_compute); + if (icompute < 0) + error->all("Compute ID for fix ave/spatial does not exist"); + if (modify->compute[icompute]->peratom_flag == 0) + error->all("Fix ave/spatial compute does not calculate per-atom info"); + nvalues = compute_size_peratom = modify->compute[icompute]->size_peratom; + if (nvalues == 0) nvalues = 1; + } + + if (me == 0) { + fprintf(fp,"Spatial-averaged data for fix %s, group %s, and %s %s\n", + id,group->names[igroup],arg[9],arg[10]); + fprintf(fp,"TimeStep Number-of-layers (one per snapshot)\n"); + fprintf(fp,"Layer Coord Atoms Value(s) (one per layer)\n"); + } + + nsum = nlayers = maxlayer = 0; + coord = NULL; + count_one = count_many = count_total = NULL; + values_one = values_many = values_total = NULL; +} + +/* ---------------------------------------------------------------------- */ + +FixAveSpatial::~FixAveSpatial() +{ + if (which == COMPUTE) delete [] id_compute; + if (me == 0) fclose(fp); + + memory->sfree(coord); + memory->sfree(count_one); + memory->sfree(count_many); + memory->sfree(count_total); + memory->destroy_2d_double_array(values_one); + memory->destroy_2d_double_array(values_many); + memory->destroy_2d_double_array(values_total); +} + +/* ---------------------------------------------------------------------- */ + +int FixAveSpatial::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveSpatial::init() +{ + // set ptrs to current compute and precompute + + if (which == COMPUTE) { + int icompute = modify->find_compute(id_compute); + if (icompute < 0) + error->all("Compute ID for fix ave/spatial does not exist"); + compute = modify->compute[icompute]; + + if (compute->id_pre) { + icompute = modify->find_compute(compute->id_pre); + if (icompute < 0) + error->all("Precompute ID for fix ave/spatial does not exist"); + precompute = modify->compute[icompute]; + } else precompute = NULL; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixAveSpatial::end_of_step() +{ + int i,j,m,ilayer,nstide; + double lo,hi; + + // if computing the first sample, setup layers + // compute current origin = boundary for some layer + // lo = layer boundary immediately below boxlo + // hi = layer boundary immediately above boxhi + // allocate and initialize arrays based on new layer count + + if (nsum == 0) { + double *boxlo = domain->boxlo; + double *boxhi = domain->boxhi; + if (originflag == LOWER) origin = boxlo[dim]; + else if (originflag == UPPER) origin = boxhi[dim]; + else if (originflag == CENTER) + origin = 0.5 * (boxlo[dim] + boxhi[dim]); + + if (origin < domain->boxlo[dim]) { + m = static_cast ((domain->boxlo[dim] - origin) * invdelta); + lo = origin + m*delta; + } else { + m = static_cast ((origin - domain->boxlo[dim]) * invdelta); + lo = origin - m*delta; + if (lo > domain->boxlo[dim]) lo -= delta; + } + if (origin < domain->boxhi[dim]) { + m = static_cast ((domain->boxhi[dim] - origin) * invdelta); + hi = origin + m*delta; + if (hi < boxhi[dim]) hi += delta; + } else { + m = static_cast ((origin - domain->boxhi[dim]) * invdelta); + hi = origin - m*delta; + } + + offset = lo; + nlayers = static_cast ((hi-lo) * invdelta + 0.5); + double volume = domain->xprd * domain->yprd * domain->zprd; + layer_volume = delta * volume/domain->prd[dim]; + + if (nlayers > maxlayer) { + maxlayer = nlayers; + coord = (double *) memory->srealloc(coord,nlayers*sizeof(double), + "ave/spatial:coord"); + count_one = (double *) + memory->srealloc(count_one,nlayers*sizeof(double), + "ave/spatial:count_one"); + count_many = (double *) + memory->srealloc(count_many,nlayers*sizeof(double), + "ave/spatial:count_many"); + count_total = (double *) + memory->srealloc(count_total,nlayers*sizeof(double), + "ave/spatial:count_total"); + values_one = memory->grow_2d_double_array(values_one,nlayers,nvalues, + "ave/spatial:values_one"); + values_many = memory->grow_2d_double_array(values_many,nlayers,nvalues, + "ave/spatial:values_many"); + values_total = memory->grow_2d_double_array(values_total,nlayers,nvalues, + "ave/spatial:values_total"); + } + + for (m = 0; m < nlayers; m++) { + coord[m] = offset + (m+0.5)*delta; + count_many[m] = count_total[m] = 0.0; + for (i = 0; i < nvalues; i++) values_many[m][i] = 0.0; + } + } + + // zero out arrays for one sample + + nsum++; + for (m = 0; m < nlayers; m++) { + count_one[m] = 0.0; + for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0; + } + + // perform the computation for one sample + // sum within each layer, only include atoms in fix group + // DENSITY_MASS adds mass + // DENSITY_NUM adds 1 to values + // ATOM adds atom vector to values + // COMPUTE adds its vector to values + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (which == DENSITY_MASS) { + int *type = atom->type; + double *mass = atom->mass; + double *rmass = atom->rmass; + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ilayer = static_cast ((x[dim][i] - offset) * invdelta); + if (ilayer < 0) ilayer = 0; + if (ilayer >= nlayers) ilayer = nlayers-1; + count_one[ilayer] += 1.0; + if (mass) values_one[ilayer][0] += mass[type[i]]; + else values_one[ilayer][0] += rmass[i]; + } + } + + } else if (which == DENSITY_NUM) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ilayer = static_cast ((x[dim][i] - offset) * invdelta); + if (ilayer < 0) ilayer = 0; + if (ilayer >= nlayers) ilayer = nlayers-1; + count_one[ilayer] += 1.0; + values_one[ilayer][0] += 1.0; + } + } + + } else if (which != COMPUTE) { + double *vector; + int nstride = 3; + if (which == VX) vector = &atom->v[0][0]; + else if (which == VY) vector = &atom->v[0][1]; + else if (which == VZ) vector = &atom->v[0][2]; + else if (which == FX) vector = &atom->f[0][0]; + else if (which == FY) vector = &atom->f[0][1]; + else if (which == FZ) vector = &atom->f[0][2]; + + m = 0; + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ilayer = static_cast ((x[dim][i] - offset) * invdelta); + if (ilayer < 0) ilayer = 0; + if (ilayer >= nlayers) ilayer = nlayers-1; + count_one[ilayer] += 1.0; + values_one[ilayer][0] += vector[m]; + } + m += nstride; + } + + } else { + if (precompute) precompute->compute_peratom(); + compute->compute_peratom(); + double *scalar = compute->scalar_atom; + double **vector = compute->vector_atom; + + m = 0; + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + ilayer = static_cast ((x[dim][i] - offset) * invdelta); + if (ilayer < 0) ilayer = 0; + if (ilayer >= nlayers) ilayer = nlayers-1; + count_one[ilayer] += 1.0; + if (compute_size_peratom == 0) values_one[ilayer][0] += scalar[i]; + else + for (j = 0; j < nvalues; j++) + values_one[ilayer][j] += vector[i][j]; + } + } + } + + // average a single sample + + if (normflag == ALL) { + for (m = 0; m < nlayers; m++) { + count_many[m] += count_one[m]; + for (j = 0; j < nvalues; j++) + values_many[m][j] += values_one[m][j]; + } + } else { + MPI_Allreduce(count_one,count_many,nlayers,MPI_DOUBLE,MPI_SUM,world); + for (m = 0; m < nlayers; m++) { + if (count_many[m] > 0.0) + for (j = 0; j < nvalues; j++) + values_many[m][j] += values_one[m][j]/count_many[m]; + count_total[m] += count_many[m]; + } + } + + // output the results + // time average across samples + // if density, also normalize by volume + + if (update->ntimestep % nfreq == 0) { + if (normflag == ALL) { + MPI_Allreduce(count_many,count_total,nlayers,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&values_many[0][0],&values_total[0][0],nlayers*nvalues, + MPI_DOUBLE,MPI_SUM,world); + for (m = 0; m < nlayers; m++) { + if (count_total[m] > 0.0) + for (j = 0; j < nvalues; j++) + values_total[m][j] /= count_total[m]; + count_total[m] /= nsum; + } + } else { + MPI_Allreduce(&values_many[0][0],&values_total[0][0],nlayers*nvalues, + MPI_DOUBLE,MPI_SUM,world); + for (m = 0; m < nlayers; m++) { + for (j = 0; j < nvalues; j++) + values_total[m][j] /= nsum; + count_total[m] /= nsum; + } + } + + if (which == DENSITY_MASS || which == DENSITY_NUM) { + for (m = 0; m < nlayers; m++) + values_total[m][0] *= count_total[m] / layer_volume; + } + + if (me == 0) { + fprintf(fp,"%d %d\n",update->ntimestep,nlayers); + for (m = 0; m < nlayers; m++) { + fprintf(fp," %d %g %g",m+1,coord[m],count_total[m]); + for (i = 0; i < nvalues; i++) fprintf(fp," %g",values_total[m][i]); + fprintf(fp,"\n"); + } + } + + nsum = 0; + } +} diff --git a/src/fix_ave_spatial.h b/src/fix_ave_spatial.h new file mode 100644 index 0000000000..bd2ff9906a --- /dev/null +++ b/src/fix_ave_spatial.h @@ -0,0 +1,52 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef FIX_AVE_SPATIAL_H +#define FIX_AVE_SPATIAL_H + +#include "stdio.h" +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAveSpatial : public Fix { + public: + FixAveSpatial(class LAMMPS *, int, char **); + ~FixAveSpatial(); + int setmask(); + void init(); + void end_of_step(); + + private: + int me; + int nfreq; + int dim,originflag,which,normflag; + double origin,delta; + char *id_compute; + FILE *fp; + + int nlayers,nvalues,nsum,maxlayer; + int compute_size_peratom; + double xscale,yscale,zscale; + double layer_volume; + double *coord; + double *count_one,*count_many,*count_total; + double **values_one,**values_many,**values_total; + double offset,invdelta; + class Compute *compute; + class Compute *precompute; +}; + +} + +#endif diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp new file mode 100644 index 0000000000..95b2ef1d69 --- /dev/null +++ b/src/fix_ave_time.cpp @@ -0,0 +1,166 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Pieter in't Veld (SNL) +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "fix_ave_time.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "group.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 8) error->all("Illegal fix ave/time command"); + + nevery = atoi(arg[3]); + nfreq = atoi(arg[4]); + + int n = strlen(arg[5]) + 1; + id_compute = new char[n]; + strcpy(id_compute,arg[5]); + + int flag = atoi(arg[6]); + + MPI_Comm_rank(world,&me); + if (me == 0) { + fp = fopen(arg[7],"w"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix ave/time file %s",arg[7]); + error->one(str); + } + } + + // setup and error check + + if (nevery <= 0) error->all("Illegal fix ave/time command"); + if (nfreq < nevery || nfreq % nevery) + error->all("Illegal fix ave/time command"); + + int icompute = modify->find_compute(id_compute); + if (icompute < 0) error->all("Compute ID for fix ave/time does not exist"); + + if (flag < 0 || flag > 2) error->all("Illegal fix ave/time command"); + sflag = vflag = 0; + if (flag == 0 || flag == 2) sflag = 1; + if (flag == 1 || flag == 2) vflag = 1; + + if (sflag && modify->compute[icompute]->scalar_flag == 0) + error->all("Fix ave/time compute does not calculate a scalar"); + if (vflag && modify->compute[icompute]->vector_flag == 0) + error->all("Fix ave/time compute does not calculate a vector"); + + if (modify->compute[icompute]->pressflag) pressure_every = nevery; + + if (me == 0) { + fprintf(fp,"Time-averaged data for fix %s, group %s, and compute %s\n", + id,group->names[modify->compute[icompute]->igroup],id_compute); + if (sflag and !vflag) + fprintf(fp,"TimeStep Value\n"); + else if (!sflag and vflag) + fprintf(fp,"TimeStep Vector-values\n"); + else if (!sflag and vflag) + fprintf(fp,"TimeStep Scalar-value Vector-values\n"); + } + + nsum = 0; + scalar = 0.0; + vector = NULL; + if (vflag) { + size_vector = modify->compute[icompute]->size_vector; + vector = new double[size_vector]; + for (int i = 0; i < size_vector; i++) vector[i] = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +FixAveTime::~FixAveTime() +{ + delete [] id_compute; + if (me == 0) fclose(fp); + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +int FixAveTime::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveTime::init() +{ + // set ptrs to current compute and precompute + + int icompute = modify->find_compute(id_compute); + if (icompute < 0) error->all("Compute ID for fix ave/time does not exist"); + compute = modify->compute[icompute]; + + if (compute->id_pre) { + icompute = modify->find_compute(compute->id_pre); + if (icompute < 0) + error->all("Precompute ID for fix ave/time does not exist"); + precompute = modify->compute[icompute]; + } else precompute = NULL; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveTime::end_of_step() +{ + int i; + + if (precompute) { + if (sflag) double tmp = precompute->compute_scalar(); + if (vflag) precompute->compute_vector(); + } + + nsum++; + if (sflag) scalar += compute->compute_scalar(); + if (vflag) { + compute->compute_vector(); + double *cvector = compute->vector; + for (i = 0; i < size_vector; i++) vector[i] += cvector[i]; + } + + if (update->ntimestep % nfreq == 0) { + if (me == 0) { + fprintf(fp,"%d",update->ntimestep); + if (sflag) fprintf(fp," %g",scalar/nsum); + if (vflag) + for (i = 0; i < size_vector; i++) fprintf(fp," %g",vector[i]/nsum); + fprintf(fp,"\n"); + } + + nsum = 0; + scalar = 0.0; + if (vflag) + for (i = 0; i < size_vector; i++) vector[i] = 0.0; + } +} diff --git a/src/fix_ave_time.h b/src/fix_ave_time.h new file mode 100644 index 0000000000..8c732064d9 --- /dev/null +++ b/src/fix_ave_time.h @@ -0,0 +1,46 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef FIX_AVE_TIME_H +#define FIX_AVE_TIME_H + +#include "stdio.h" +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAveTime : public Fix { + public: + FixAveTime(class LAMMPS *, int, char **); + ~FixAveTime(); + int setmask(); + void init(); + void end_of_step(); + + private: + int me; + + int nfreq; + char *id_compute; + FILE *fp; + + int sflag,vflag; + int size_vector,nsum; + double scalar,*vector; + class Compute *compute; + class Compute *precompute; +}; + +} + +#endif diff --git a/src/min_cg.cpp b/src/min_cg.cpp index f21a98825b..7e2dc2dd7c 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -72,12 +72,12 @@ void MinCG::init() delete [] fixarg; fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1]; - // zero gradient vectors before first atom exchange + // zero local vectors before first atom exchange - setup_vectors(); + set_local_vectors(); for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0; - // virial_thermo is how virial should be computed on thermo timesteps + // virial_thermo = how virial computed on thermo timesteps // 1 = computed explicity by pair, 2 = computed implicitly by pair if (force->newton_pair) virial_thermo = 2; @@ -127,28 +127,31 @@ void MinCG::init() void MinCG::run() { - double tmp,*f; + int i; + double tmp; // set initial force & energy setup(); - setup_vectors(); output->thermo->compute_pe(); - ecurrent = output->thermo->potential_energy; + energy = output->thermo->potential_energy; + energy += energy_extra; // stats for Finish to print - einitial = ecurrent; + einitial = energy; - f = atom->f[0]; tmp = 0.0; - for (int i = 0; i < ndof; i++) tmp += f[i]*f[i]; + for (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 (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp); + for (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 @@ -169,16 +172,19 @@ void MinCG::run() output->next_dump_any = update->ntimestep; if (output->restart_every) output->next_restart = update->ntimestep; output->next_thermo = update->ntimestep; - int ntmp; - double *xtmp,*htmp,etmp; - eng_force(&ntmp,&xtmp,&htmp,&etmp); + eng_force(); 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 @@ -188,17 +194,20 @@ void MinCG::run() // stats for Finish to print - efinal = ecurrent; + efinal = energy; f = atom->f[0]; tmp = 0.0; - for (int i = 0; i < ndof; i++) tmp += f[i]*f[i]; + for (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 (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp); + for (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])); } /* ---------------------------------------------------------------------- @@ -207,6 +216,20 @@ 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 @@ -222,7 +245,7 @@ void MinCG::setup() comm->borders(); neighbor->build(); neighbor->ncalls = 0; - setup_vectors(); + set_local_vectors(); // compute all forces @@ -247,6 +270,7 @@ void MinCG::setup() if (force->newton) comm->reverse_communicate(); modify->setup(); + energy_extra = modify->min_energy(xextra,fextra); output->setup(1); } @@ -259,14 +283,14 @@ void MinCG::iterate(int n) { int i,gradsearch,fail; double alpha,beta,gg,dot[2],dotall[2]; - double *f; - f = atom->f[0]; - for (int i = 0; i < ndof; i++) h[i] = g[i] = f[i]; + for (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.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; @@ -277,8 +301,8 @@ void MinCG::iterate(int n) // line minimization along direction h from current atom->x - eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval); + eprevious = energy; + fail = (this->*linemin)(neval); // if max_eval exceeded, all done // if linemin failed or energy did not decrease sufficiently: @@ -287,8 +311,8 @@ void MinCG::iterate(int n) if (neval >= update->max_eval) break; - if (fail || fabs(ecurrent-eprevious) <= - update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) { + if (fail || fabs(energy-eprevious) <= + update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS)) { if (gradsearch == 1) break; gradsearch = -1; } @@ -299,14 +323,17 @@ 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; @@ -319,6 +346,10 @@ 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 @@ -334,9 +365,11 @@ void MinCG::iterate(int n) set ndof and vector pointers after atoms have migrated ------------------------------------------------------------------------- */ -void MinCG::setup_vectors() +void MinCG::set_local_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]; @@ -346,11 +379,11 @@ void MinCG::setup_vectors() /* ---------------------------------------------------------------------- evaluate potential energy and forces may migrate atoms - new energy stored in ecurrent and returned (in case caller not in class) - negative gradient will be stored in atom->f + energy = new objective function energy = poteng of atoms + eng_extra + atom->f, fextra = negative gradient of objective function ------------------------------------------------------------------------- */ -void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng) +void MinCG::eng_force() { // check for reneighboring // always communicate since minimizer moved atoms @@ -375,7 +408,7 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng) timer->stamp(TIME_COMM); neighbor->build(); timer->stamp(TIME_NEIGHBOR); - setup_vectors(); + set_local_vectors(); } // eflag is always set, since minimizer needs potential energy @@ -409,21 +442,17 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng) timer->stamp(TIME_COMM); } - // fixes that affect minimization + // min_post_force = forces on atoms 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(); - ecurrent = output->thermo->potential_energy; - - // return updated ptrs to caller since atoms may have migrated - - *pndof = ndof; - *px = atom->x[0]; - *ph = h; - *peng = ecurrent; + energy = output->thermo->potential_energy; + energy += energy_extra; } /* ---------------------------------------------------------------------- @@ -487,25 +516,15 @@ void MinCG::force_clear(int vflag) /* ---------------------------------------------------------------------- line minimization methods - find minimum-energy starting at x along dir direction - input: n = # of degrees of freedom on this proc - x = ptr to atom->x[0] as vector - dir = search direction as vector - eng = current energy at initial x - min/max dist = min/max distance to move any atom coord - output: return 0 if successful move, set alpha - return 1 if failed, no move, no need to set alpha - 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 + find minimum-energy starting at x along h direction + update atom->x by alpha, call eng_force() for result + alpha = distance moved along h to set x to minimun-energy configuration + return 0 if successful move, 1 if failed (no move) + insure last call to eng_force() is consistent with return + if fail, eng_force() of original x + if succeed, eng_force() at x + alpha*h + eng_force() may migrate atoms due to neighbor list build + therefore linemin routines CANNOT store atom-based quantities ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- @@ -515,52 +534,55 @@ void MinCG::force_clear(int vflag) quit as soon as energy starts to rise ------------------------------------------------------------------------- */ -int MinCG::linemin_scan(int n, double *x, double *dir, double eng, - double mindist, double maxdist, - double &alpha, int &nfunc) +int MinCG::linemin_scan(int &nfunc) { int i; - double fmax,fme,elowest,alphamin,alphamax,alphalast; + double fmax,fme,elowest,alpha,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 < n; i++) fme = MAX(fme,fabs(dir[i])); + for (i = 0; i < ndof; i++) fme = MAX(fme,fabs(h[i])); MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world); - if (fmax == 0.0) return 1; + for (i = 0; i < ndof_extra; i++) fmax = MAX(fmax,fabs(hextra[i])); - alphamin = mindist/fmax; - alphamax = maxdist/fmax; + if (fmax == 0.0) return 1; + alphamin = dmin/fmax; + alphamax = dmax/fmax; // if minstep is already uphill, fail // if eng increases, stop and return previous alpha // if alphamax, stop and return alphamax - elowest = eng; + elowest = energy; alpha = alphamin; while (1) { - for (i = 0; i < n; i++) x[i] += alpha*dir[i]; - eng_force(&n,&x,&dir,&eng); + for (i = 0; i < ndof; i++) x[i] += alpha*h[i]; + for (i = 0; i < ndof_extra; i++) xextra[i] += alpha*hextra[i]; + eng_force(); nfunc++; - if (alpha == alphamin && eng >= elowest) { - for (i = 0; i < n; i++) x[i] -= alpha*dir[i]; - eng_force(&n,&x,&dir,&eng); + if (alpha == alphamin && energy >= elowest) { + for (i = 0; i < ndof; i++) x[i] -= alpha*h[i]; + for (i = 0; i < ndof_extra; i++) xextra[i] -= alpha*hextra[i]; + eng_force(); nfunc++; return 1; } - if (eng > elowest) { - for (i = 0; i < n; i++) x[i] += (alphalast-alpha)*dir[i]; - eng_force(&n,&x,&dir,&eng); + if (energy > elowest) { + for (i = 0; i < ndof; i++) x[i] += (alphalast-alpha)*h[i]; + for (i = 0; i < ndof_extra; i++) + xextra[i] += (alphalast-alpha)*hextra[i]; + eng_force(); nfunc++; alpha = alphalast; return 0; } if (alpha == alphamax) return 0; - elowest = eng; + elowest = energy; alphalast = alpha; alpha *= SCAN_FACTOR; if (alpha > alphamax) alpha = alphamax; @@ -574,43 +596,44 @@ int MinCG::linemin_scan(int n, double *x, double *dir, double eng, prevents successvive func evals further apart in x than maxdist ------------------------------------------------------------------------- */ -int MinCG::linemin_secant(int n, double *x, double *dir, double eng, - double mindist, double maxdist, - double &alpha, int &nfunc) +int MinCG::linemin_secant(int &nfunc) { int i,iter; - double eta,eta_prev,sigma0,sigmamax,alphadelta,fme,fmax,dsq,e0,tmp; - double *f; + double eta,eta_prev,sigma0,sigmamax,alpha,alphadelta,fme,fmax,dsq,e0,tmp; double epssq = SECANT_EPS * SECANT_EPS; // stopping criterion for secant iterations fme = 0.0; - for (i = 0; i < n; i++) fme += dir[i]*dir[i]; + for (i = 0; i < ndof; i++) fme += h[i]*h[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 < n; i++) fme = MAX(fme,fabs(dir[i])); + for (i = 0; i < ndof; i++) fme = MAX(fme,fabs(h[i])); MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world); - if (fmax == 0.0) return 1; + for (i = 0; i < ndof_extra; i++) fmax = MAX(fmax,fabs(hextra[i])); - sigma0 = mindist/fmax; - sigmamax = maxdist/fmax; + if (fmax == 0.0) return 1; + sigma0 = dmin/fmax; + sigmamax = dmax/fmax; // eval func at sigma0 // test if minstep is already uphill - e0 = eng; - for (i = 0; i < n; i++) x[i] += sigma0*dir[i]; - eng_force(&n,&x,&dir,&eng); + e0 = energy; + for (i = 0; i < ndof; i++) x[i] += sigma0*h[i]; + for (i = 0; i < ndof_extra; i++) xextra[i] += sigma0*hextra[i]; + eng_force(); nfunc++; - if (eng >= e0) { - for (i = 0; i < n; i++) x[i] -= sigma0*dir[i]; - eng_force(&n,&x,&dir,&eng); + if (energy >= e0) { + for (i = 0; i < ndof; i++) x[i] -= sigma0*h[i]; + for (i = 0; i < ndof_extra; i++) xextra[i] -= sigma0*hextra[i]; + eng_force(); nfunc++; return 1; } @@ -618,24 +641,25 @@ int MinCG::linemin_secant(int n, double *x, double *dir, double eng, // secant iterations // alphadelta = new increment to move, alpha = accumulated move - f = atom->f[0]; tmp = 0.0; - for (i = 0; i < n; i++) tmp -= f[i]*dir[i]; + for (i = 0; i < ndof; i++) tmp -= f[i]*h[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 < n; i++) x[i] += alphadelta*dir[i]; - eng_force(&n,&x,&dir,&eng); + for (i = 0; i < ndof; i++) x[i] += alphadelta*h[i]; + for (i = 0; i < ndof_extra; i++) xextra[i] += alphadelta*hextra[i]; + eng_force(); nfunc++; - f = atom->f[0]; tmp = 0.0; - for (i = 0; i < n; i++) tmp -= f[i]*dir[i]; + for (i = 0; i < ndof; i++) tmp -= f[i]*h[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; diff --git a/src/min_cg.h b/src/min_cg.h index f8fcfaaf9e..4db5c15378 100644 --- a/src/min_cg.h +++ b/src/min_cg.h @@ -27,32 +27,39 @@ class MinCG : public Min { virtual void iterate(int); protected: - int virial_thermo; // what vflag should be on thermo steps (1,2) - int pairflag,torqueflag,granflag; + int virial_thermo; // what vflag should be on thermo steps (1,2) + int pairflag,torqueflag,granflag; // force clear flags 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 ecurrent; // current potential energy - double mindist,maxdist; // min/max dist for coord delta in line search + double mindist,maxdist; // min/max dist for coord delta in line search - int ndof; // # of degrees-of-freedom on this proc - double *g,*h; // local portion of gradient, searchdir vectors + int ndof; // 3N degrees-of-freedom on this proc + double *x; // vec of 3N coords, ptr to atom->x[0] + 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 - typedef int (MinCG::*FnPtr)(int, double *, double *, double, - double, double, double &, int &); - FnPtr linemin; // ptr to linemin functions + int ndof_extra; // extra degrees of freedom to include in min + double energy_extra; // extra potential energy + double *xextra; // extra vecs associated with ndof_extra + double *fextra; + double *gextra; + double *hextra; - int linemin_scan(int, double *, double *, double, - double, double, double &, int &); - int linemin_secant(int, double *, double *, double, - double, double, double &, int &); + double energy; // potential energy of atoms and extra dof + + typedef int (MinCG::*FnPtr)(int &); + FnPtr linemin; // ptr to linemin functions + int linemin_scan(int &); + int linemin_secant(int &); void setup(); - void setup_vectors(); - void eng_force(int *, double **, double **, double *); + void set_local_vectors(); + void eng_force(); void force_clear(int); }; diff --git a/src/min_cg_fr.cpp b/src/min_cg_fr.cpp index fbac308558..666449445f 100644 --- a/src/min_cg_fr.cpp +++ b/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 = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval); + eprevious = energy; + fail = (this->*linemin)(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(ecurrent-eprevious) <= - update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) { + if (fail || fabs(energy-eprevious) <= + update->tolerance * 0.5*(fabs(energy) + 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,6 +96,10 @@ 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 diff --git a/src/min_sd.cpp b/src/min_sd.cpp index 83a678af10..36debfd64b 100644 --- a/src/min_sd.cpp +++ b/src/min_sd.cpp @@ -35,10 +35,9 @@ 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; @@ -48,28 +47,29 @@ void MinSD::iterate(int n) // line minimization along direction h from current atom->x - eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval); + eprevious = energy; + fail = (this->*linemin)(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(ecurrent-eprevious) <= - update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) + if (fail || fabs(energy-eprevious) <= + update->tolerance * 0.5*(fabs(energy) + 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 diff --git a/src/style_class2.h b/src/style_class2.h index 3ba9b32df4..e69de29bb2 100644 --- a/src/style_class2.h +++ b/src/style_class2.h @@ -1,56 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_dpd.h b/src/style_dpd.h index 8ce617c0c2..e69de29bb2 100644 --- a/src/style_dpd.h +++ b/src/style_dpd.h @@ -1,28 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_granular.h b/src/style_granular.h index 7b0f4b6a71..e69de29bb2 100644 --- a/src/style_granular.h +++ b/src/style_granular.h @@ -1,50 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_kspace.h b/src/style_kspace.h index 52cae87bdb..e69de29bb2 100644 --- a/src/style_kspace.h +++ b/src/style_kspace.h @@ -1,38 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_manybody.h b/src/style_manybody.h index 449895a530..e69de29bb2 100644 --- a/src/style_manybody.h +++ b/src/style_manybody.h @@ -1,28 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_molecule.h b/src/style_molecule.h index 6e424035a9..e69de29bb2 100644 --- a/src/style_molecule.h +++ b/src/style_molecule.h @@ -1,116 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_opt.h b/src/style_opt.h index 061816e136..e69de29bb2 100644 --- a/src/style_opt.h +++ b/src/style_opt.h @@ -1,30 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 diff --git a/src/style_xtc.h b/src/style_xtc.h index 7110dda312..e69de29bb2 100644 --- a/src/style_xtc.h +++ b/src/style_xtc.h @@ -1,20 +0,0 @@ -/* ---------------------------------------------------------------------- - 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