From ead9fe70299ab223ff88121133d4ef9ccf4ebeca Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 6 Mar 2013 16:44:51 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9582 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-OMP/respa_omp.cpp | 358 +++++++++++++++++++++++++++++++++++++ src/USER-OMP/respa_omp.h | 111 ++++++++++++ 2 files changed, 469 insertions(+) create mode 100644 src/USER-OMP/respa_omp.cpp create mode 100644 src/USER-OMP/respa_omp.h diff --git a/src/USER-OMP/respa_omp.cpp b/src/USER-OMP/respa_omp.cpp new file mode 100644 index 0000000000..b1ee76d0b0 --- /dev/null +++ b/src/USER-OMP/respa_omp.cpp @@ -0,0 +1,358 @@ +/* ---------------------------------------------------------------------- + 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 authors: Mark Stevens (SNL), Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "respa_omp.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "atom.h" +#include "force.h" +#include "pair.h" +#include "bond.h" +#include "angle.h" +#include "dihedral.h" +#include "improper.h" +#include "kspace.h" +#include "output.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "fix_respa.h" +#include "timer.h" +#include "memory.h" +#include "error.h" + +#if defined(_OPENMP) +#include +#endif + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +RespaOMP::RespaOMP(LAMMPS *lmp, int narg, char **arg) + : Respa(lmp, narg, arg),ThrOMP(lmp, THR_INTGR) +{ +} + +/* ---------------------------------------------------------------------- + initialization before run +------------------------------------------------------------------------- */ + +void RespaOMP::init() +{ + Respa::init(); + + if (atom->torque) + error->all(FLERR,"Extended particles are not supported by respa/omp\n"); +} + +/* ---------------------------------------------------------------------- + setup before run +------------------------------------------------------------------------- */ + +void RespaOMP::setup() +{ + if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n"); + + update->setupflag = 1; + + // setup domain, communication and neighboring + // acquire ghosts + // build neighbor lists + + atom->setup(); + modify->setup_pre_exchange(); + if (triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + domain->reset_box(); + comm->setup(); + if (neighbor->style) neighbor->setup_bins(); + comm->exchange(); + if (atom->sortfreq > 0) atom->sort(); + comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + domain->image_check(); + domain->box_too_small_check(); + modify->setup_pre_neighbor(); + neighbor->build(); + neighbor->ncalls = 0; + + // compute all forces + + ev_set(update->ntimestep); + + for (int ilevel = 0; ilevel < nlevels; ilevel++) { + force_clear(newton[ilevel]); + modify->setup_pre_force_respa(vflag,ilevel); + if (level_bond == ilevel && force->bond) + force->bond->compute(eflag,vflag); + if (level_angle == ilevel && force->angle) + force->angle->compute(eflag,vflag); + if (level_dihedral == ilevel && force->dihedral) + force->dihedral->compute(eflag,vflag); + if (level_improper == ilevel && force->improper) + force->improper->compute(eflag,vflag); + if (level_pair == ilevel && pair_compute_flag) + force->pair->compute(eflag,vflag); + if (level_inner == ilevel && pair_compute_flag) + force->pair->compute_inner(); + if (level_middle == ilevel && pair_compute_flag) + force->pair->compute_middle(); + if (level_outer == ilevel && pair_compute_flag) + force->pair->compute_outer(eflag,vflag); + if (level_kspace == ilevel && force->kspace) { + force->kspace->setup(); + if (kspace_compute_flag) force->kspace->compute(eflag,vflag); + } + + // reduce forces from per-thread arrays, if needed + if (!fix->get_reduced()) { + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; +#if defined(_OPENMP) +#pragma omp parallel default(none) +#endif + { +#if defined(_OPENMP) + int tid = omp_get_thread_num(); +#else + int tid = 0; +#endif + data_reduce_thr(atom->f[0], nall, nthreads, 3, tid); + } + fix->did_reduce(); + } + + if (newton[ilevel]) comm->reverse_comm(); + copy_f_flevel(ilevel); + } + + modify->setup(vflag); + sum_flevel_f(); + output->setup(); + update->setupflag = 0; +} + +/* ---------------------------------------------------------------------- + setup without output + flag = 0 = just force calculation + flag = 1 = reneighbor and force calculation +------------------------------------------------------------------------- */ + +void RespaOMP::setup_minimal(int flag) +{ + update->setupflag = 1; + + // setup domain, communication and neighboring + // acquire ghosts + // build neighbor lists + + if (flag) { + modify->setup_pre_exchange(); + if (triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + domain->reset_box(); + comm->setup(); + if (neighbor->style) neighbor->setup_bins(); + comm->exchange(); + comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + domain->image_check(); + domain->box_too_small_check(); + modify->setup_pre_neighbor(); + neighbor->build(); + neighbor->ncalls = 0; + } + + // compute all forces + + ev_set(update->ntimestep); + + for (int ilevel = 0; ilevel < nlevels; ilevel++) { + force_clear(newton[ilevel]); + modify->setup_pre_force_respa(vflag,ilevel); + if (level_bond == ilevel && force->bond) + force->bond->compute(eflag,vflag); + if (level_angle == ilevel && force->angle) + force->angle->compute(eflag,vflag); + if (level_dihedral == ilevel && force->dihedral) + force->dihedral->compute(eflag,vflag); + if (level_improper == ilevel && force->improper) + force->improper->compute(eflag,vflag); + if (level_pair == ilevel && pair_compute_flag) + force->pair->compute(eflag,vflag); + if (level_inner == ilevel && pair_compute_flag) + force->pair->compute_inner(); + if (level_middle == ilevel && pair_compute_flag) + force->pair->compute_middle(); + if (level_outer == ilevel && pair_compute_flag) + force->pair->compute_outer(eflag,vflag); + if (level_kspace == ilevel && force->kspace) { + force->kspace->setup(); + if (kspace_compute_flag) force->kspace->compute(eflag,vflag); + } + + // reduce forces from per-thread arrays, if needed + if (!fix->get_reduced()) { + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; +#if defined(_OPENMP) +#pragma omp parallel default(none) +#endif + { +#if defined(_OPENMP) + int tid = omp_get_thread_num(); +#else + int tid = 0; +#endif + data_reduce_thr(atom->f[0], nall, nthreads, 3, tid); + } + fix->did_reduce(); + } + + if (newton[ilevel]) comm->reverse_comm(); + copy_f_flevel(ilevel); + } + + modify->setup(vflag); + sum_flevel_f(); + update->setupflag = 0; +} + +/* ---------------------------------------------------------------------- */ + +void RespaOMP::recurse(int ilevel) +{ + copy_flevel_f(ilevel); + + for (int iloop = 0; iloop < loop[ilevel]; iloop++) { + + modify->initial_integrate_respa(vflag,ilevel,iloop); + if (modify->n_post_integrate_respa) + modify->post_integrate_respa(ilevel,iloop); + + if (ilevel) recurse(ilevel-1); + + // at outermost level, check on rebuilding neighbor list + // at innermost level, communicate + // at middle levels, do nothing + + if (ilevel == nlevels-1) { + int nflag = neighbor->decide(); + if (nflag) { + if (modify->n_pre_exchange) modify->pre_exchange(); + if (triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + if (domain->box_change) { + domain->reset_box(); + comm->setup(); + if (neighbor->style) neighbor->setup_bins(); + } + timer->stamp(); + comm->exchange(); + if (atom->sortfreq > 0 && + update->ntimestep >= atom->nextsort) atom->sort(); + comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + timer->stamp(TIME_COMM); + if (modify->n_pre_neighbor) modify->pre_neighbor(); + neighbor->build(); + timer->stamp(TIME_NEIGHBOR); + } + + } else if (ilevel == 0) { + timer->stamp(); + comm->forward_comm(); + timer->stamp(TIME_COMM); + } + + force_clear(newton[ilevel]); + if (modify->n_pre_force_respa) + modify->pre_force_respa(vflag,ilevel,iloop); + + timer->stamp(); + if (level_bond == ilevel && force->bond) { + force->bond->compute(eflag,vflag); + timer->stamp(TIME_BOND); + } + if (level_angle == ilevel && force->angle) { + force->angle->compute(eflag,vflag); + timer->stamp(TIME_BOND); + } + if (level_dihedral == ilevel && force->dihedral) { + force->dihedral->compute(eflag,vflag); + timer->stamp(TIME_BOND); + } + if (level_improper == ilevel && force->improper) { + force->improper->compute(eflag,vflag); + timer->stamp(TIME_BOND); + } + if (level_pair == ilevel && pair_compute_flag) { + force->pair->compute(eflag,vflag); + timer->stamp(TIME_PAIR); + } + if (level_inner == ilevel && pair_compute_flag) { + force->pair->compute_inner(); + timer->stamp(TIME_PAIR); + } + if (level_middle == ilevel && pair_compute_flag) { + force->pair->compute_middle(); + timer->stamp(TIME_PAIR); + } + if (level_outer == ilevel && pair_compute_flag) { + force->pair->compute_outer(eflag,vflag); + timer->stamp(TIME_PAIR); + } + if (level_kspace == ilevel && kspace_compute_flag) { + force->kspace->compute(eflag,vflag); + timer->stamp(TIME_KSPACE); + } + + // reduce forces from per-thread arrays, if needed + if (!fix->get_reduced()) { + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; +#if defined(_OPENMP) +#pragma omp parallel default(none) +#endif + { +#if defined(_OPENMP) + int tid = omp_get_thread_num(); +#else + int tid = 0; +#endif + data_reduce_thr(atom->f[0], nall, nthreads, 3, tid); + } + fix->did_reduce(); + } + + if (newton[ilevel]) { + comm->reverse_comm(); + timer->stamp(TIME_COMM); + } + + if (modify->n_post_force_respa) + modify->post_force_respa(vflag,ilevel,iloop); + modify->final_integrate_respa(ilevel,iloop); + } + + copy_f_flevel(ilevel); +} + diff --git a/src/USER-OMP/respa_omp.h b/src/USER-OMP/respa_omp.h new file mode 100644 index 0000000000..89ddbe78af --- /dev/null +++ b/src/USER-OMP/respa_omp.h @@ -0,0 +1,111 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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 INTEGRATE_CLASS + +IntegrateStyle(respa/omp,RespaOMP) + +#else + +#ifndef LMP_RESPA_OMP_H +#define LMP_RESPA_OMP_H + +#include "respa.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class RespaOMP : public Respa, public ThrOMP { + public: + RespaOMP(class LAMMPS *, int, char **); + virtual ~RespaOMP() {} + virtual void init(); + virtual void setup(); + virtual void setup_minimal(int); + + protected: + virtual void recurse(int); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Respa levels must be >= 1 + +Self-explanatory. + +E: Cannot set both respa pair and inner/middle/outer + +In the rRESPA integrator, you must compute pairwise potentials either +all together (pair), or in pieces (inner/middle/outer). You can't do +both. + +E: Must set both respa inner and outer + +Cannot use just the inner or outer option with respa without using the +other. + +E: Cannot set respa middle without inner/outer + +In the rRESPA integrator, you must define both a inner and outer +setting in order to use a middle setting. + +E: Invalid order of forces within respa levels + +For respa, ordering of force computations within respa levels must +obey certain rules. E.g. bonds cannot be compute less frequently than +angles, pairwise forces cannot be computed less frequently than +kspace, etc. + +W: One or more respa levels compute no forces + +This is computationally inefficient. + +E: Respa inner cutoffs are invalid + +The first cutoff must be <= the second cutoff. + +E: Respa middle cutoffs are invalid + +The first cutoff must be <= the second cutoff. + +W: No fixes defined, atoms won't move + +If you are not using a fix like nve, nvt, npt then atom velocities and +coordinates will not be updated during timestepping. + +W: Fix shake with rRESPA computes invalid pressures + +This is a known bug in LAMMPS that has not yet been fixed. If you use +SHAKE with rRESPA and perform a constant volume simulation (e.g. using +fix npt) this only affects the output pressure, not the dynamics of +the simulation. If you use SHAKE with rRESPA and perform a constant +pressure simulation (e.g. using fix npt) then you will be +equilibrating to the wrong volume. + +E: Pair style does not support rRESPA inner/middle/outer + +You are attempting to use rRESPA options with a pair style that +does not support them. + +*/