diff --git a/doc/src/Eqs/pair_ufm.jpg b/doc/src/Eqs/pair_ufm.jpg new file mode 100644 index 0000000000..40273da680 Binary files /dev/null and b/doc/src/Eqs/pair_ufm.jpg differ diff --git a/doc/src/Eqs/pair_ufm.tex b/doc/src/Eqs/pair_ufm.tex new file mode 100644 index 0000000000..8a6dde5520 --- /dev/null +++ b/doc/src/Eqs/pair_ufm.tex @@ -0,0 +1,14 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ + E = -\varepsilon\, \ln{\left[1-\exp{\left(-r^{2}/\sigma^{2}\right)}\right]} \qquad r < r_c +$$ + +$$ + \varepsilon = p\,k_B\,T +$$ + +\end{document} + diff --git a/doc/src/pair_ufm.txt b/doc/src/pair_ufm.txt new file mode 100644 index 0000000000..2be35b0d4b --- /dev/null +++ b/doc/src/pair_ufm.txt @@ -0,0 +1,135 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +pair_style ufm command :h3 +pair_style ufm/gpu command :h3 +pair_style ufm/omp command :h3 +pair_style ufm/opt command :h3 + +[Syntax:] + +pair_style ufm cutoff :pre + +cutoff = global cutoff for {ufm} interactions (distance units) :ul + +[Examples:] + +pair_style ufm 4.0 +pair_coeff 1 1 100.0 1.0 2.5 +pair_coeff * * 100.0 1.0 :pre + + +pair_style ufm 4.0 +pair_coeff * * 10.0 1.0 +variable prefactor equal ramp(10,100) +fix 1 all adapt 1 pair ufm epsilon * * v_prefactor :pre + +[Description:] + +Style {ufm} computes pairwise interactions using the Uhlenbeck-Ford model (UFM) potential "(Paula Leite2016)"_#PL2 which is given by + +:c,image(Eqs/pair_ufm.jpg) + +where rc is the cutoff, sigma is a distance-scale and epsilon is an energy-scale, i.e., a product of Boltzmann constant kB, temperature T and the Uhlenbeck-Ford p-parameter which is responsible +to control the softness of the interactions "(Paula Leite2017)"_#PL1. +This model is useful as a reference system for fluid-phase free-energy calculations "(Paula Leite2016)"_#PL2. + +The following coefficients must be defined for each pair of atom types +via the "pair_coeff"_pair_coeff.html command as in the examples above, +or in the data file or restart files read by the +"read_data"_read_data.html or "read_restart"_read_restart.html +commands, or by mixing as described below: + +epsilon (energy units) +sigma (distance units) +cutoff (distance units) :ul + +The last coefficient is optional. If not specified, the global {ufm} +cutoff is used. + + +The "fix adapt"_fix_adapt.html command can be used to vary epsilon and sigma for this pair style over the course of a simulation, in which case +pair_coeff settings for epsilon and sigma must still be specified, but will be +overridden. For example these commands will vary the prefactor epsilon for +all pairwise interactions from 10.0 at the beginning to 100.0 at the end +of a run: + +variable prefactor equal ramp(10,100) +fix 1 all adapt 1 pair ufm epsilon * * v_prefactor :pre + +NOTE: The thermodynamic integration procedure can be performed with this potential using "fix adapt"_fix_adapt.html. This command will rescale the force on each atom by varying a scale variable, which always starts with value 1.0. The syntax is the same described above, however, changing epsilon to scale. A detailed explanation of how to use this command and perform nonequilibrium thermodynamic integration in LAMMPS is given in the paper by "(Freitas)"_#Freitas. + +:line + +Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are +functionally the same as the corresponding style without the suffix. +They have been optimized to run faster, depending on your available +hardware, as discussed in "Section 5"_Section_accelerate.html +of the manual. The accelerated styles take the same arguments and +should produce the same results, except for round-off and precision +issues. + +These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, +USER-OMP and OPT packages, respectively. They are only enabled if +LAMMPS was built with those packages. See the "Making +LAMMPS"_Section_start.html#start_3 section for more info. + +You can specify the accelerated styles explicitly in your input script +by including their suffix, or you can use the "-suffix command-line +switch"_Section_start.html#start_6 when you invoke LAMMPS, or you can +use the "suffix"_suffix.html command in your input script. + +See "Section 5"_Section_accelerate.html of the manual for +more instructions on how to use the accelerated styles effectively. + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +For atom type pairs I,J and I != J, the A coefficient and cutoff +distance for this pair style can be mixed. A is always mixed via a +{geometric} rule. The cutoff is mixed according to the pair_modify +mix value. The default mix value is {geometric}. See the +"pair_modify" command for details. + +This pair style support the "pair_modify"_pair_modify.html shift option for the energy of the pair interaction. + +The "pair_modify"_pair_modify.html table and tail are not relevant for this +pair style. + +This pair style does not support the "pair_modify"_pair_modify.html tail option for adding long-range tail corrections to energy and pressure. + +This pair style writes its information to "binary restart +files"_restart.html, so pair_style and pair_coeff commands do not need +to be specified in an input script that reads a restart file. + +This pair style can only be used via the {pair} keyword of the +"run_style respa"_run_style.html command. It does not support the +{inner}, {middle}, {outer} keywords. + +:line + +[Restrictions:] none + +[Related commands:] + +"pair_coeff"_pair_coeff.html, "fix adapt"_fix_adapt.html + +[Default:] none + + +:link(PL1) +[(Paula Leite2017)] Paula Leite, Santos-Florez, and de Koning, Phys Rev E, 96, +32115 (2017). + +:link(PL2) +[(Paula Leite2016)] Paula Leite , Freitas, Azevedo, and de Koning, J Chem Phys, 126, +044509 (2016). + +:link(Freitas) +[(Freitas)] Freitas, Asta, and de Koning, Computational Materials Science, 112, 333 (2016). \ No newline at end of file diff --git a/lib/gpu/Nvidia.makefile b/lib/gpu/Nvidia.makefile index ee2eb72632..5f692cf66c 100644 --- a/lib/gpu/Nvidia.makefile +++ b/lib/gpu/Nvidia.makefile @@ -76,7 +76,8 @@ OBJS = $(OBJ_DIR)/lal_atom.o $(OBJ_DIR)/lal_ans.o \ $(OBJ_DIR)/lal_coul.o $(OBJ_DIR)/lal_coul_ext.o \ $(OBJ_DIR)/lal_coul_debye.o $(OBJ_DIR)/lal_coul_debye_ext.o \ $(OBJ_DIR)/lal_zbl.o $(OBJ_DIR)/lal_zbl_ext.o \ - $(OBJ_DIR)/lal_lj_cubic.o $(OBJ_DIR)/lal_lj_cubic_ext.o + $(OBJ_DIR)/lal_lj_cubic.o $(OBJ_DIR)/lal_lj_cubic_ext.o \ + $(OBJ_DIR)/lal_ufm.o $(OBJ_DIR)/lal_ufm_ext.o CBNS = $(OBJ_DIR)/device.cubin $(OBJ_DIR)/device_cubin.h \ $(OBJ_DIR)/atom.cubin $(OBJ_DIR)/atom_cubin.h \ @@ -131,7 +132,8 @@ CBNS = $(OBJ_DIR)/device.cubin $(OBJ_DIR)/device_cubin.h \ $(OBJ_DIR)/coul.cubin $(OBJ_DIR)/coul_cubin.h \ $(OBJ_DIR)/coul_debye.cubin $(OBJ_DIR)/coul_debye_cubin.h \ $(OBJ_DIR)/zbl.cubin $(OBJ_DIR)/zbl_cubin.h \ - $(OBJ_DIR)/lj_cubic.cubin $(OBJ_DIR)/lj_cubic_cubin.h + $(OBJ_DIR)/lj_cubic.cubin $(OBJ_DIR)/lj_cubic_cubin.h \ + $(OBJ_DIR)/ufm.cubin $(OBJ_DIR)/ufm_cubin.h all: $(OBJ_DIR) $(GPU_LIB) $(EXECS) @@ -705,6 +707,18 @@ $(OBJ_DIR)/dpd.cubin: lal_dpd.cu lal_precision.h lal_preprocessor.h $(OBJ_DIR)/dpd_cubin.h: $(OBJ_DIR)/dpd.cubin $(OBJ_DIR)/dpd.cubin $(BIN2C) -c -n dpd $(OBJ_DIR)/dpd.cubin > $(OBJ_DIR)/dpd_cubin.h +$(OBJ_DIR)/ufm.cubin: lal_ufm.cu lal_precision.h lal_preprocessor.h + $(CUDA) --cubin -DNV_KERNEL -o $@ lal_ufm.cu + +$(OBJ_DIR)/ufm_cubin.h: $(OBJ_DIR)/ufm.cubin $(OBJ_DIR)/ufm.cubin + $(BIN2C) -c -n ufm $(OBJ_DIR)/ufm.cubin > $(OBJ_DIR)/ufm_cubin.h + +$(OBJ_DIR)/lal_ufm.o: $(ALL_H) lal_ufm.h lal_ufm.cpp $(OBJ_DIR)/ufm_cubin.h $(OBJ_DIR)/lal_base_atomic.o + $(CUDR) -o $@ -c lal_ufm.cpp -I$(OBJ_DIR) + +$(OBJ_DIR)/lal_ufm_ext.o: $(ALL_H) lal_ufm.h lal_ufm_ext.cpp lal_base_atomic.h + $(CUDR) -o $@ -c lal_ufm_ext.cpp -I$(OBJ_DIR) + $(OBJ_DIR)/lal_dpd.o: $(ALL_H) lal_dpd.h lal_dpd.cpp $(OBJ_DIR)/dpd_cubin.h $(OBJ_DIR)/lal_base_dpd.o $(CUDR) -o $@ -c lal_dpd.cpp -I$(OBJ_DIR) diff --git a/lib/gpu/README b/lib/gpu/README index b26897e885..15b65516ac 100644 --- a/lib/gpu/README +++ b/lib/gpu/README @@ -135,6 +135,7 @@ Current styles supporting GPU acceleration: 38 yukawa/colloid 39 yukawa 40 pppm + 41 ufm MULTIPLE LAMMPS PROCESSES diff --git a/lib/gpu/lal_ufm.cpp b/lib/gpu/lal_ufm.cpp new file mode 100644 index 0000000000..c7aa2cca39 --- /dev/null +++ b/lib/gpu/lal_ufm.cpp @@ -0,0 +1,172 @@ +/*************************************************************************** + ufm.cpp + ------------------- + Rodolfo Paula Leite (Unicamp/Brazil) + Maurice de Koning (Unicamp/Brazil) + + Class for acceleration of the ufm pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : pl.rodolfo@gmail.com + dekoning@ifi.unicamp.br + ***************************************************************************/ + +#if defined(USE_OPENCL) +#include "ufm_cl.h" +#elif defined(USE_CUDART) +const char *ufm=0; +#else +#include "ufm_cubin.h" +#endif + +#include "lal_ufm.h" +#include +using namespace LAMMPS_AL; +#define UFMT UFM + +extern Device device; + +template +UFMT::UFM() : BaseAtomic(), _allocated(false) { +} + +template +UFMT::~UFM() { + clear(); +} + +template +int UFMT::bytes_per_atom(const int max_nbors) const { + return this->bytes_per_atom_atomic(max_nbors); +} + +template +int UFMT::init(const int ntypes, + double **host_cutsq, double **host_uf1, + double **host_uf2, double **host_uf3, + double **host_uf4, double **host_offset, + double *host_special_lj, const int nlocal, + const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *_screen) { + int success; + success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split, + _screen,ufm,"k_ufm"); + if (success!=0) + return success; + + // If atom type constants fit in shared memory use fast kernel + int lj_types=ntypes; + shared_types=false; + int max_shared_types=this->device->max_shared_types(); + if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) { + lj_types=max_shared_types; + shared_types=true; + } + _lj_types=lj_types; + + // Allocate a host write buffer for data initialization + UCL_H_Vec host_write(lj_types*lj_types*32,*(this->ucl_device), + UCL_WRITE_ONLY); + + for (int i=0; iucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,uf1,host_write,host_uf1,host_uf2, + host_cutsq); + + uf3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY); + this->atom->type_pack4(ntypes,lj_types,uf3,host_write,host_uf3,host_uf4, + host_offset); + + UCL_H_Vec dview; + sp_lj.alloc(4,*(this->ucl_device),UCL_READ_ONLY); + dview.view(host_special_lj,4,*(this->ucl_device)); + ucl_copy(sp_lj,dview,false); + + _allocated=true; + this->_max_bytes=uf1.row_bytes()+uf3.row_bytes()+sp_lj.row_bytes(); + return 0; +} + +template +void UFMT::reinit(const int ntypes, double **host_cutsq, double **host_uf1, + double **host_uf2, double **host_uf3, + double **host_uf4, double **host_offset) { + // Allocate a host write buffer for data initialization + UCL_H_Vec host_write(_lj_types*_lj_types*32,*(this->ucl_device), + UCL_WRITE_ONLY); + + for (int i=0; i<_lj_types*_lj_types; i++) + host_write[i]=0.0; + + this->atom->type_pack4(ntypes,_lj_types,uf1,host_write,host_uf1,host_uf2, + host_cutsq); + this->atom->type_pack4(ntypes,_lj_types,uf3,host_write,host_uf3,host_uf4, + host_offset); +} + +template +void UFMT::clear() { + if (!_allocated) + return; + _allocated=false; + + uf1.clear(); + uf3.clear(); + sp_lj.clear(); + this->clear_atomic(); +} + +template +double UFMT::host_memory_usage() const { + return this->host_memory_usage_atomic()+sizeof(UFM); +} + +// --------------------------------------------------------------------------- +// Calculate energies, forces, and torques +// --------------------------------------------------------------------------- +template +void UFMT::loop(const bool _eflag, const bool _vflag) { + // Compute the block size and grid size to keep all cores busy + const int BX=this->block_size(); + int eflag, vflag; + if (_eflag) + eflag=1; + else + eflag=0; + + if (_vflag) + vflag=1; + else + vflag=0; + + int GX=static_cast(ceil(static_cast(this->ans->inum())/ + (BX/this->_threads_per_atom))); + + int ainum=this->ans->inum(); + int nbor_pitch=this->nbor->nbor_pitch(); + this->time_pair.start(); + if (shared_types) { + this->k_pair_fast.set_size(GX,BX); + this->k_pair_fast.run(&this->atom->x, &uf1, &uf3, &sp_lj, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, + &vflag, &ainum, &nbor_pitch, + &this->_threads_per_atom); + } else { + this->k_pair.set_size(GX,BX); + this->k_pair.run(&this->atom->x, &uf1, &uf3, &_lj_types, &sp_lj, + &this->nbor->dev_nbor, &this->_nbor_data->begin(), + &this->ans->force, &this->ans->engv, &eflag, &vflag, + &ainum, &nbor_pitch, &this->_threads_per_atom); + } + this->time_pair.stop(); +} + +template class UFM; diff --git a/lib/gpu/lal_ufm.cu b/lib/gpu/lal_ufm.cu new file mode 100644 index 0000000000..51c4df3b5b --- /dev/null +++ b/lib/gpu/lal_ufm.cu @@ -0,0 +1,188 @@ +/*************************************************************************** + ufm.cu + ------------------- + Rodolfo Paula Leite (Unicamp/Brazil) + Maurice de Koning (Unicamp/Brazil) + + Device code for acceleration of the ufm pair style + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : pl.rodolfo@gmail.com + dekoning@ifi.unicamp.br + ***************************************************************************/ + +#ifdef NV_KERNEL +#include "lal_aux_fun1.h" +#ifndef _DOUBLE_DOUBLE +texture pos_tex; +#else +texture pos_tex; +#endif +#else +#define pos_tex x_ +#endif + +__kernel void k_ufm(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict uf1, + const __global numtyp4 *restrict uf3, + const int lj_types, + const __global numtyp *restrict sp_lj, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, const int inum, + const int nbor_pitch, const int t_per_atom) { + int tid, ii, offset; + atom_info(t_per_atom,ii,tid,offset); + + acctyp energy=(acctyp)0; + acctyp4 f; + f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0; + acctyp virial[6]; + for (int i=0; i<6; i++) + virial[i]=(acctyp)0; + + if (ii0) { + energy += - factor_lj * uf3[mtype].x*log(1.0 - expuf) - uf3[mtype].z; + } + if (vflag>0) { + virial[0] += delx*delx*force; + virial[1] += dely*dely*force; + virial[2] += delz*delz*force; + virial[3] += delx*dely*force; + virial[4] += delx*delz*force; + virial[5] += dely*delz*force; + } + } + + } // for nbor + store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, + ans,engv); + } // if ii +} + +__kernel void k_ufm_fast(const __global numtyp4 *restrict x_, + const __global numtyp4 *restrict uf1_in, + const __global numtyp4 *restrict uf3_in, + const __global numtyp *restrict sp_lj_in, + const __global int * dev_nbor, + const __global int * dev_packed, + __global acctyp4 *restrict ans, + __global acctyp *restrict engv, + const int eflag, const int vflag, const int inum, + const int nbor_pitch, const int t_per_atom) { + int tid, ii, offset; + atom_info(t_per_atom,ii,tid,offset); + + __local numtyp4 uf1[MAX_SHARED_TYPES*MAX_SHARED_TYPES]; + __local numtyp4 uf3[MAX_SHARED_TYPES*MAX_SHARED_TYPES]; + __local numtyp sp_lj[4]; + if (tid<4) + sp_lj[tid]=sp_lj_in[tid]; + if (tid0) + uf3[tid]=uf3_in[tid]; + } + + acctyp energy=(acctyp)0; + acctyp4 f; + f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0; + acctyp virial[6]; + for (int i=0; i<6; i++) + virial[i]=(acctyp)0; + + __syncthreads(); + + if (ii0) { + energy += - factor_lj * uf3[mtype].x * log(1.0 - expuf) - uf3[mtype].z; + } + if (vflag>0) { + virial[0] += delx*delx*force; + virial[1] += dely*dely*force; + virial[2] += delz*delz*force; + virial[3] += delx*dely*force; + virial[4] += delx*delz*force; + virial[5] += dely*delz*force; + } + } + + } // for nbor + store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, + ans,engv); + } // if ii +} + diff --git a/lib/gpu/lal_ufm.h b/lib/gpu/lal_ufm.h new file mode 100644 index 0000000000..aeeaacbe99 --- /dev/null +++ b/lib/gpu/lal_ufm.h @@ -0,0 +1,86 @@ +/*************************************************************************** + ufm.h + ------------------- + Rodolfo Paula Leite (Unicamp/Brazil) + Maurice de Koning (Unicamp/Brazil) + + Class for acceleration of the ufm pair style. + + __________________________________________________________________________ + This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) + __________________________________________________________________________ + + begin : + email : pl.rodolfo@gmail.com + dekoning@ifi.unicamp.br + ***************************************************************************/ + +#ifndef LAL_UFM_H +#define LAL_UFM_H + +#include "lal_base_atomic.h" + +namespace LAMMPS_AL { + +template +class UFM : public BaseAtomic { + public: + UFM(); + ~UFM(); + + /// Clear any previous data and set up for a new LAMMPS run + /** \param max_nbors initial number of rows in the neighbor matrix + * \param cell_size cutoff + skin + * \param gpu_split fraction of particles handled by device + * + * Returns: + * - 0 if successfull + * - -1 if fix gpu not found + * - -3 if there is an out of memory error + * - -4 if the GPU library was not compiled for GPU + * - -5 Double precision is not supported on card **/ + int init(const int ntypes, double **host_cutsq, + double **host_uf1, double **host_uf2, double **host_uf3, + double **host_uf4, double **host_offset, double *host_special_lj, + const int nlocal, const int nall, const int max_nbors, + const int maxspecial, const double cell_size, + const double gpu_split, FILE *screen); + + /// Send updated coeffs from host to device (to be compatible with fix adapt) + void reinit(const int ntypes, double **host_cutsq, + double **host_uf1, double **host_uf2, double **host_uf3, + double **host_uf4, double **host_offset); + + /// Clear all host and device data + /** \note This is called at the beginning of the init() routine **/ + void clear(); + + /// Returns memory usage on device per atom + int bytes_per_atom(const int max_nbors) const; + + /// Total host memory used by library for pair style + double host_memory_usage() const; + + // --------------------------- TYPE DATA -------------------------- + + /// uf1.x = uf1, uf1.y = uf2, uf1.z = cutsq + UCL_D_Vec uf1; + /// uf3.x = uf3, uf3.y = uf4, uf3.z = offset + UCL_D_Vec uf3; + /// Special LJ values + UCL_D_Vec sp_lj; + + /// If atom type constants fit in shared memory, use fast kernels + bool shared_types; + + /// Number of atom types + int _lj_types; + + private: + bool _allocated; + void loop(const bool _eflag, const bool _vflag); +}; + +} + +#endif diff --git a/src/GPU/Install.sh b/src/GPU/Install.sh index f4aeaa2706..88f47a3dc4 100644 --- a/src/GPU/Install.sh +++ b/src/GPU/Install.sh @@ -131,6 +131,8 @@ action pair_zbl_gpu.cpp action pair_zbl_gpu.h action pppm_gpu.cpp pppm.cpp action pppm_gpu.h pppm.cpp +action pair_ufm_gpu.cpp +action pair_ufm_gpu.h # edit 2 Makefile.package files to include/exclude package info diff --git a/src/GPU/pair_ufm_gpu.cpp b/src/GPU/pair_ufm_gpu.cpp new file mode 100644 index 0000000000..96af0dc069 --- /dev/null +++ b/src/GPU/pair_ufm_gpu.cpp @@ -0,0 +1,240 @@ +/* -*- 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. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#include +#include +#include +#include "pair_ufm_gpu.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "integrate.h" +#include "memory.h" +#include "error.h" +#include "neigh_request.h" +#include "universe.h" +#include "update.h" +#include "domain.h" +#include +#include "gpu_extra.h" + +using namespace LAMMPS_NS; + +// External functions from cuda library for atom decomposition + +int ufml_gpu_init(const int ntypes, double **cutsq, double **host_uf1, + double **host_uf2, double **host_uf3, double **host_uf4, + double **offset, double *special_lj, const int nlocal, + const int nall, const int max_nbors, const int maxspecial, + const double cell_size, int &gpu_mode, FILE *screen); + +int ufml_gpu_reinit(const int ntypes, double **cutsq, double **host_uf1, + double **host_uf2, double **host_uf3, double **host_uf4, + double **offset); + +void ufml_gpu_clear(); +int ** ufml_gpu_compute_n(const int ago, const int inum, + const int nall, double **host_x, int *host_type, + double *sublo, double *subhi, tagint *tag, int **nspecial, + tagint **special, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + int **ilist, int **jnum, + const double cpu_time, bool &success); +void ufml_gpu_compute(const int ago, const int inum, const int nall, + double **host_x, int *host_type, int *ilist, int *numj, + int **firstneigh, const bool eflag, const bool vflag, + const bool eatom, const bool vatom, int &host_start, + const double cpu_time, bool &success); +double ufml_gpu_bytes(); + +/* ---------------------------------------------------------------------- */ + +PairUFMGPU::PairUFMGPU(LAMMPS *lmp) : PairUFM(lmp), gpu_mode(GPU_FORCE) +{ + respa_enable = 0; + cpu_time = 0.0; + GPU_EXTRA::gpu_ready(lmp->modify, lmp->error); +} + +/* ---------------------------------------------------------------------- + free all arrays +------------------------------------------------------------------------- */ + +PairUFMGPU::~PairUFMGPU() +{ + ufml_gpu_clear(); +} + +/* ---------------------------------------------------------------------- */ + +void PairUFMGPU::compute(int eflag, int vflag) +{ + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + int nall = atom->nlocal + atom->nghost; + int inum, host_start; + + bool success = true; + int *ilist, *numneigh, **firstneigh; + if (gpu_mode != GPU_FORCE) { + inum = atom->nlocal; + firstneigh = ufml_gpu_compute_n(neighbor->ago, inum, nall, + atom->x, atom->type, domain->sublo, + domain->subhi, atom->tag, atom->nspecial, + atom->special, eflag, vflag, eflag_atom, + vflag_atom, host_start, + &ilist, &numneigh, cpu_time, success); + } else { + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + ufml_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type, + ilist, numneigh, firstneigh, eflag, vflag, eflag_atom, + vflag_atom, host_start, cpu_time, success); + } + if (!success) + error->one(FLERR,"Insufficient memory on accelerator"); + + if (host_startnewton_pair) + error->all(FLERR,"Cannot use newton pair with ufm/gpu pair style"); + + // Repeat cutsq calculation because done after call to init_style + double maxcut = -1.0; + double cut; + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = i; j <= atom->ntypes; j++) { + if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) { + cut = init_one(i,j); + cut *= cut; + if (cut > maxcut) + maxcut = cut; + cutsq[i][j] = cutsq[j][i] = cut; + } else + cutsq[i][j] = cutsq[j][i] = 0.0; + } + } + double cell_size = sqrt(maxcut) + neighbor->skin; + + int maxspecial=0; + if (atom->molecular) + maxspecial=atom->maxspecial; + int success = ufml_gpu_init(atom->ntypes+1, cutsq, uf1, uf2, uf3, uf4, + offset, force->special_lj, atom->nlocal, + atom->nlocal+atom->nghost, 300, maxspecial, + cell_size, gpu_mode, screen); + GPU_EXTRA::check_flag(success,error,world); + + if (gpu_mode == GPU_FORCE) { + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairUFMGPU::reinit() +{ + Pair::reinit(); + + ufml_gpu_reinit(atom->ntypes+1, cutsq, uf1, uf2, uf3, uf4, offset); +} + +/* ---------------------------------------------------------------------- */ + +double PairUFMGPU::memory_usage() +{ + double bytes = Pair::memory_usage(); + return bytes + ufml_gpu_bytes(); +} + +/* ---------------------------------------------------------------------- */ + +void PairUFMGPU::cpu_compute(int start, int inum, int eflag, int vflag, + int *ilist, int *numneigh, int **firstneigh) { + int i,j,ii,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,expuf,factor_lj; + int *jlist; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + double *special_lj = force->special_lj; + + + // loop over neighbors of my atoms + + for (ii = start; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + expuf = exp(- rsq * uf2[itype][jtype]); + fpair = factor_lj * uf1[itype][jtype] * expuf /(1.0 - expuf); + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + + if (eflag) { + evdwl = -factor_lj * uf3[itype][jtype] * log(1.0 - expuf) - offset[itype][jtype]; + } + + if (evflag) ev_tally_full(i,evdwl,0.0,fpair,delx,dely,delz); + } + } + } +} diff --git a/src/GPU/pair_ufm_gpu.h b/src/GPU/pair_ufm_gpu.h new file mode 100644 index 0000000000..59b883f3aa --- /dev/null +++ b/src/GPU/pair_ufm_gpu.h @@ -0,0 +1,65 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(ufm/gpu,PairUFMGPU) + +#else + +#ifndef LMP_PAIR_UFM_GPU_H +#define LMP_PAIR_UFM_GPU_H + +#include "pair_ufm.h" + +namespace LAMMPS_NS { + +class PairUFMGPU : public PairUFM { + public: + PairUFMGPU(LAMMPS *lmp); + ~PairUFMGPU(); + void cpu_compute(int, int, int, int, int *, int *, int **); + void compute(int, int); + void init_style(); + void reinit(); + double memory_usage(); + + enum { GPU_FORCE, GPU_NEIGH, GPU_HYB_NEIGH }; + + private: + int gpu_mode; + double cpu_time; +}; + +} +#endif +#endif + +/* ERROR/WARNING messages: + +E: Insufficient memory on accelerator + +There is insufficient memory on one of the devices specified for the gpu +package + +E: Cannot use newton pair with ufm/gpu pair style + +Self-explanatory. + +*/ diff --git a/src/OPT/Install.sh b/src/OPT/Install.sh index ca1231c615..c6ae2b914b 100644 --- a/src/OPT/Install.sh +++ b/src/OPT/Install.sh @@ -46,3 +46,5 @@ action pair_lj_long_coul_long_opt.cpp pair_lj_long_coul_long.cpp action pair_lj_long_coul_long_opt.h pair_lj_long_coul_long.cpp action pair_morse_opt.cpp action pair_morse_opt.h +action pair_ufm_opt.cpp +action pair_ufm_opt.h diff --git a/src/OPT/pair_ufm_opt.cpp b/src/OPT/pair_ufm_opt.cpp new file mode 100644 index 0000000000..1cf504674d --- /dev/null +++ b/src/OPT/pair_ufm_opt.cpp @@ -0,0 +1,197 @@ +/* ---------------------------------------------------------------------- + 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: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#include +#include "pair_ufm_opt.h" +#include "atom.h" +#include "force.h" +#include "neigh_list.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairUFMOpt::PairUFMOpt(LAMMPS *lmp) : PairUFM(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void PairUFMOpt::compute(int eflag, int vflag) +{ + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + if (evflag) { + if (eflag) { + if (force->newton_pair) return eval<1,1,1>(); + else return eval<1,1,0>(); + } else { + if (force->newton_pair) return eval<1,0,1>(); + else return eval<1,0,0>(); + } + } else { + if (force->newton_pair) return eval<0,0,1>(); + else return eval<0,0,0>(); + } +} + +/* ---------------------------------------------------------------------- */ + +template < int EVFLAG, int EFLAG, int NEWTON_PAIR > +void PairUFMOpt::eval() +{ + typedef struct { double x,y,z; } vec3_t; + + typedef struct { + double cutsq,uf1,uf2,uf3,uf4,scale,offset; + double _pad[2]; + } fast_alpha_t; + + int i,j,ii,jj,inum,jnum,itype,jtype,sbindex; + double factor_lj; + double evdwl = 0.0; + + double** _noalias x = atom->x; + double** _noalias f = atom->f; + int* _noalias type = atom->type; + int nlocal = atom->nlocal; + double* _noalias special_lj = force->special_lj; + + inum = list->inum; + int* _noalias ilist = list->ilist; + int** _noalias firstneigh = list->firstneigh; + int* _noalias numneigh = list->numneigh; + + vec3_t* _noalias xx = (vec3_t*)x[0]; + vec3_t* _noalias ff = (vec3_t*)f[0]; + + int ntypes = atom->ntypes; + int ntypes2 = ntypes*ntypes; + + fast_alpha_t* _noalias fast_alpha = + (fast_alpha_t*) malloc(ntypes2*sizeof(fast_alpha_t)); + for (i = 0; i < ntypes; i++) for (j = 0; j < ntypes; j++) { + fast_alpha_t& a = fast_alpha[i*ntypes+j]; + a.cutsq = cutsq[i+1][j+1]; + a.uf1 = uf1[i+1][j+1]; + a.uf2 = uf2[i+1][j+1]; + a.uf3 = uf3[i+1][j+1]; + a.uf4 = uf4[i+1][j+1]; + a.scale = scale[i+1][j+1]; + a.offset = offset[i+1][j+1]; + } + fast_alpha_t* _noalias tabsix = fast_alpha; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + double xtmp = xx[i].x; + double ytmp = xx[i].y; + double ztmp = xx[i].z; + itype = type[i] - 1; + int* _noalias jlist = firstneigh[i]; + jnum = numneigh[i]; + + double tmpfx = 0.0; + double tmpfy = 0.0; + double tmpfz = 0.0; + + fast_alpha_t* _noalias tabsixi = (fast_alpha_t*)&tabsix[itype*ntypes]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + sbindex = sbmask(j); + + if (sbindex == 0) { + double delx = xtmp - xx[j].x; + double dely = ytmp - xx[j].y; + double delz = ztmp - xx[j].z; + double rsq = delx*delx + dely*dely + delz*delz; + + jtype = type[j] - 1; + + fast_alpha_t& a = tabsixi[jtype]; + + if (rsq < a.cutsq) { + double expuf = exp(- rsq * a.uf2); + double fpair = a.scale * a.uf1 * expuf / (1.0 - expuf); + + tmpfx += delx*fpair; + tmpfy += dely*fpair; + tmpfz += delz*fpair; + if (NEWTON_PAIR || j < nlocal) { + ff[j].x -= delx*fpair; + ff[j].y -= dely*fpair; + ff[j].z -= delz*fpair; + } + + if (EFLAG) evdwl = - a.uf3 * log(1.0 - expuf) - a.offset; + + if (EVFLAG) + ev_tally(i,j,nlocal,NEWTON_PAIR, + evdwl,0.0,fpair,delx,dely,delz); + } + + } else { + factor_lj = special_lj[sbindex]; + j &= NEIGHMASK; + + double delx = xtmp - xx[j].x; + double dely = ytmp - xx[j].y; + double delz = ztmp - xx[j].z; + double rsq = delx*delx + dely*dely + delz*delz; + + int jtype1 = type[j]; + jtype = jtype1 - 1; + + fast_alpha_t& a = tabsixi[jtype]; + if (rsq < a.cutsq) { + fast_alpha_t& a = tabsixi[jtype]; + double expuf = exp(- rsq * a.uf2); + double fpair = a.scale * factor_lj * a.uf1 * expuf / (1.0 - expuf); + + tmpfx += delx*fpair; + tmpfy += dely*fpair; + tmpfz += delz*fpair; + if (NEWTON_PAIR || j < nlocal) { + ff[j].x -= delx*fpair; + ff[j].y -= dely*fpair; + ff[j].z -= delz*fpair; + } + + if (EFLAG) { + evdwl = - a.uf3 * log(1.0 - expuf) - a.offset; + evdwl *= factor_lj; + } + + if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + ff[i].x += tmpfx; + ff[i].y += tmpfy; + ff[i].z += tmpfz; + } + + free(fast_alpha); fast_alpha = 0; + + if (vflag_fdotr) virial_fdotr_compute(); +} diff --git a/src/OPT/pair_ufm_opt.h b/src/OPT/pair_ufm_opt.h new file mode 100644 index 0000000000..edac708403 --- /dev/null +++ b/src/OPT/pair_ufm_opt.h @@ -0,0 +1,45 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(ufm/opt,PairUFMOpt) + +#else + +#ifndef LMP_PAIR_UFM_OPT_H +#define LMP_PAIR_UFM_OPT_H + +#include "pair_ufm.h" + +namespace LAMMPS_NS { + +class PairUFMOpt : public PairUFM { + public: + PairUFMOpt(class LAMMPS *); + void compute(int, int); + + private: + template < int EVFLAG, int EFLAG, int NEWTON_PAIR > void eval(); +}; + +} + +#endif +#endif diff --git a/src/USER-OMP/pair_ufm_omp.cpp b/src/USER-OMP/pair_ufm_omp.cpp new file mode 100644 index 0000000000..b2e2cd29ee --- /dev/null +++ b/src/USER-OMP/pair_ufm_omp.cpp @@ -0,0 +1,159 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#include +#include "pair_ufm_omp.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" + +#include "suffix.h" +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairUFMOMP::PairUFMOMP(LAMMPS *lmp) : + PairUFM(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; + respa_enable = 0; +} + +/* ---------------------------------------------------------------------- */ + +void PairUFMOMP::compute(int eflag, int vflag) +{ + if (eflag || vflag) { + ev_setup(eflag,vflag); + } else evflag = vflag_fdotr = 0; + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(eflag,vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); + + if (evflag) { + if (eflag) { + if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); + else eval<1,1,0>(ifrom, ito, thr); + } else { + if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); + else eval<1,0,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); + else eval<0,0,0>(ifrom, ito, thr); + } + + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +template +void PairUFMOMP::eval(int iifrom, int iito, ThrData * const thr) +{ + int i,j,ii,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,expuf,factor; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + + const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; + dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const int * _noalias const type = atom->type; + const int nlocal = atom->nlocal; + const double * _noalias const special_lj = force->special_lj; + double fxtmp,fytmp,fztmp; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = iifrom; ii < iito; ++ii) { + + i = ilist[ii]; + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + fxtmp=fytmp=fztmp=0.0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + expuf = exp(- rsq * uf2[itype][jtype]); + fpair = factor * scale[itype][jtype] * uf1[itype][jtype] * expuf /(1.0 - expuf); + + fxtmp += delx*fpair; + fytmp += dely*fpair; + fztmp += delz*fpair; + if (NEWTON_PAIR || j < nlocal) { + f[j].x -= delx*fpair; + f[j].y -= dely*fpair; + f[j].z -= delz*fpair; + } + + if (EFLAG) { + evdwl = -uf3[itype][jtype] * log(1.0 - expuf) - offset[itype][jtype]; + evdwl *= factor; + } + + if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR, + evdwl,0.0,fpair,delx,dely,delz,thr); + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairUFMOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairUFM::memory_usage(); + + return bytes; +} diff --git a/src/USER-OMP/pair_ufm_omp.h b/src/USER-OMP/pair_ufm_omp.h new file mode 100644 index 0000000000..2a01da15d0 --- /dev/null +++ b/src/USER-OMP/pair_ufm_omp.h @@ -0,0 +1,50 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(ufm/omp,PairUFMOMP) + +#else + +#ifndef LMP_PAIR_UFM_OMP_H +#define LMP_PAIR_UFM_OMP_H + +#include "pair_ufm.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairUFMOMP : public PairUFM, public ThrOMP { + + public: + PairUFMOMP(class LAMMPS *); + + virtual void compute(int, int); + virtual double memory_usage(); + + private: + template + void eval(int ifrom, int ito, ThrData * const thr); +}; + +} + +#endif +#endif diff --git a/src/pair_ufm.cpp b/src/pair_ufm.cpp new file mode 100644 index 0000000000..5307ced365 --- /dev/null +++ b/src/pair_ufm.cpp @@ -0,0 +1,376 @@ +/* -*- 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. + ------------------------------------------------------------------------- */ + +/* ----------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_ufm.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairUFM::PairUFM(LAMMPS *lmp) : Pair(lmp) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairUFM::~PairUFM() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cut); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(scale); + memory->destroy(uf1); + memory->destroy(uf2); + memory->destroy(uf3); + memory->destroy(uf4); + memory->destroy(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairUFM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq, expuf, factor; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + expuf = exp(- rsq * uf2[itype][jtype]); + fpair = factor * scale[itype][jtype] * uf1[itype][jtype] * expuf /(1.0 - expuf); + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + evdwl = -uf3[itype][jtype] * log(1.0 - expuf) - offset[itype][jtype]; + evdwl *= factor; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairUFM::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(scale,n+1,n+1,"pair:scale"); + memory->create(uf1,n+1,n+1,"pair:uf1"); + memory->create(uf2,n+1,n+1,"pair:uf2"); + memory->create(uf3,n+1,n+1,"pair:uf3"); + memory->create(uf4,n+1,n+1,"pair:uf4"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairUFM::settings(int narg, char **arg) +{ + if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + + cut_global = force->numeric(FLERR,arg[0]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairUFM::coeff(int narg, char **arg) +{ + if (narg < 4 || narg > 5) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + + double cut_one = cut_global; + + if (narg == 5) cut_one = force->numeric(FLERR,arg[4]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + epsilon[i][j] = epsilon_one; + sigma[i][j] = sigma_one; + scale[i][j] = 1.0; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairUFM::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + uf1[i][j] = 2.0 * epsilon[i][j] / pow(sigma[i][j],2.0); + uf2[i][j] = 1.0 / pow(sigma[i][j],2.0); + uf3[i][j] = epsilon[i][j]; + uf4[i][j] = sigma[i][j]; + + if (offset_flag) { + double ratio = pow(cut[i][j] / sigma[i][j],2.0); + offset[i][j] = - epsilon[i][j] * log ( 1.0 - exp( -ratio )) ; + } else offset[i][j] = 0.0; + + uf1[j][i] = uf1[i][j]; + uf2[j][i] = uf2[i][j]; + uf3[j][i] = uf3[i][j]; + uf4[j][i] = uf4[i][j]; + scale[j][i] = scale[i][j]; + offset[j][i] = offset[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairUFM::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&epsilon[i][j],sizeof(double),1,fp); + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairUFM::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&epsilon[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairUFM::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairUFM::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairUFM::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairUFM::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairUFM::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double expuf,phiuf; + expuf = exp(- rsq * uf2[itype][jtype]); + fforce = factor_lj * uf1[itype][jtype] * expuf /(1.0 - expuf); + phiuf = - uf3[itype][jtype] * log(1.0 - expuf) - offset[itype][jtype]; + return factor_lj * phiuf; +} + +/* ---------------------------------------------------------------------- */ + +void *PairUFM::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"scale") == 0) return (void *) scale; + return NULL; +} diff --git a/src/pair_ufm.h b/src/pair_ufm.h new file mode 100644 index 0000000000..2161c2acaf --- /dev/null +++ b/src/pair_ufm.h @@ -0,0 +1,76 @@ +/* -*- 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. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: + Rodolfo Paula Leite (Unicamp/Brazil) - pl.rodolfo@gmail.com + Maurice de Koning (Unicamp/Brazil) - dekoning@ifi.unicamp.br + ------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(ufm,PairUFM) + +#else + +#ifndef LMP_PAIR_UFM_H +#define LMP_PAIR_UFM_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairUFM : public Pair { + public: + PairUFM(class LAMMPS *); + virtual ~PairUFM(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void write_data(FILE *); + void write_data_all(FILE *); + double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double cut_global; + double **cut,**scale; + double **epsilon,**sigma; + double **uf1,**uf2,**uf3,**uf4,**offset; + + virtual void allocate(); +}; + +} + +#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: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +*/