/* ---------------------------------------------------------------------- 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 "ctype.h" #include "stdlib.h" #include "pair_hybrid_weighted.h" #include "atom.h" #include "force.h" #include "neighbor.h" #include "neigh_request.h" #include "memory.h" #include "modify.h" #include "update.h" #include "integrate.h" #include "compute.h" #include "error.h" using namespace LAMMPS_NS; //#define MIN(a,b) ((a) < (b) ? (a) : (b)) //#define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ PairHybridWeighted::PairHybridWeighted(LAMMPS *lmp) : PairHybridOverlay(lmp) { forces = NULL; //oldtags = NULL; nmax = 0; } PairHybridWeighted::~PairHybridWeighted() { if (allocated) { delete[] hweights; memory->destroy(forces); } } /* ---------------------------------------------------------------------- perform substyle computes, and weigh energies and forces * store forces before substyle is called * force = oldforce*(1-weight)+newforce*weight ------------------------------------------------------------------------- */ void PairHybridWeighted::compute(int eflag, int vflag) { int i,j,m,n; n = atom->nlocal; if (force->newton_pair) n += atom->nghost; // grow forces array if necessary if( n > nmax ) { memory->destroy(forces); //memory->destroy(oldtags); nmax = atom->nmax; memory->create(forces,3*nmax,"pair:oldforces"); //memory->create(oldtags,nmax,"pair:oldtags"); } // // if no_virial_compute is set and global component of incoming vflag = 2 // reset vflag as if global component were 1 // necessary since one or more sub-styles cannot compute virial as F dot r if (no_virial_fdotr_compute && vflag % 4 == 2) vflag = 1 + vflag/4 * 4; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = eflag_global = vflag_global = eflag_atom = vflag_atom = 0; // check if global component of incoming vflag = 2 // if so, reset vflag passed to substyle as if it were 0 // necessary so substyle will not invoke virial_compute() int vflag_substyle; if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4; else vflag_substyle = vflag; double **f = atom->f; class Compute *scomp = modify->compute[id_pe]; // clear final force buffer and energy memset(forces,0,3*n*sizeof(double)); eng_vdwl = 0.0; eng_coul = 0.0; int nall; if (force->newton) nall = atom->nlocal + atom->nghost; else nall = atom->nlocal; size_t nbytes = sizeof(double) * nall; for (m = 0; m < nstyles; m++) { // get weight double l = hweights[m]; // clear forces from previous substyle if( m>0 && nbytes ) { //update->integrate->force_clear(); memset(&(atom->f[0][0]),0,3*nbytes); } // compute forces of substyle styles[m]->compute(eflag,vflag_substyle); // add weighted forces if( l>0.0 ) { for( i = 0; i < n; ++i ) { // this may shuffle atoms around!!!! forces[3*i+0] += l*f[i][0]; forces[3*i+1] += l*f[i][1]; forces[3*i+2] += l*f[i][2]; //oldtags[i] = atom->tag[i]; } } // store the total energy of each style scomp->vector[m] = styles[m]->eng_vdwl + styles[m]->eng_coul; fprintf(screen,"m=%d scomp=%f (%f+%f), l=%f\n", m, styles[m]->eng_vdwl + styles[m]->eng_coul, styles[m]->eng_vdwl, styles[m]->eng_coul, l ); if (eflag_global) { eng_vdwl += styles[m]->eng_vdwl*l; eng_coul += styles[m]->eng_coul*l; } if (vflag_global) { for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n]*l; } if (eflag_atom) { double *eatom_substyle = styles[m]->eatom; for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i]*l; } if (vflag_atom) { double **vatom_substyle = styles[m]->vatom; for (i = 0; i < n; i++) for (j = 0; j < 6; j++) vatom[i][j] += vatom_substyle[i][j]*l; } } // apply final aggregated force for( i = 0; i < n; ++i ) { f[i][0] = forces[3*i+0]; f[i][1] = forces[3*i+1]; f[i][2] = forces[3*i+2]; } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- combine sub-style neigh list requests and create new ones if needed copied from hybrid overlay ------------------------------------------------------------------------- */ void PairHybridWeighted::modify_requests() { int i,j; NeighRequest *irq,*jrq; // loop over pair requests only // if a previous list is same kind with same skip attributes // then make this one a copy list of that one // works whether both lists are no-skip or yes-skip // will not point a list at a copy list, but at copy list's parent for (i = 0; i < neighbor->nrequest; i++) { if (!neighbor->requests[i]->pair) continue; irq = neighbor->requests[i]; for (j = 0; j < i; j++) { if (!neighbor->requests[j]->pair) continue; jrq = neighbor->requests[j]; if (irq->same_kind(jrq) && irq->same_skip(jrq)) { irq->copy = 1; irq->otherlist = j; break; } } } // perform same operations on skip lists as pair style = hybrid PairHybrid::modify_requests(); } // wrap the original settings routine and allocate space for the weights void PairHybridWeighted::settings(int narg, char **arg) { PairHybrid::settings(narg,arg); nhweights = nstyles; hweights = new double[nstyles]; for( int i = 0; i < nstyles; ++i ) hweights[i]=1.0; fprintf(screen,"Allocating space for %d weights.\n", nstyles); // create a compute to store the potential energies of all substyles // id=hype, compute group = all char **newarg = new char*[4]; newarg[0] = (char *) "hype"; newarg[1] = (char *) "all"; newarg[2] = (char *) "store/pe"; newarg[3] = new char[30]; snprintf( newarg[3], 30, "%d", nstyles ); modify->add_compute(4,newarg); id_pe = modify->find_compute("hype"); delete [] newarg; }