diff --git a/README b/README index 5b13bb4dec..780cffc301 100644 --- a/README +++ b/README @@ -26,6 +26,7 @@ and directories: README this file LICENSE the GNU General Public License (GPL) bench benchmark problems +couple code coupling examples, calling LAMMPS as a library doc documentation examples simple test problems lib libraries LAMMPS can be linked with diff --git a/couple/README b/couple/README new file mode 100644 index 0000000000..675ef2f795 --- /dev/null +++ b/couple/README @@ -0,0 +1,199 @@ +This directory has examples of how to use LAMMPS as a library, either +by itself or in tandem with another code or library. + +This directory is meant to illustrate what is possible when coupling +codes or calling LAMMPS as a library. All of the examples provided +are just for demonstration purposes. The physics they calculate is +too simple for modeling a realistic problem. + +See these sections of the LAMMPS manaul for details: + +2.4 Building LAMMPS as a library (doc/Section_start.html#2_4) +4.10 Coupling LAMMPS to other codes (doc/Section_howto.html#4_10) + +In all of the examples included here, LAMMPS must first be built as a +library. Basically, you type something like + +make makelib +make -f Makefile.lib g++ + +in the LAMMPS src directory to create liblmp_g++.a + +The library interface to LAMMPS is in src/library.cpp. Routines can +be easily added to this file so an external program can perform the +LAMMPS tasks desired. + +-------------------------------- + +These are the sub-directories included in this directory: + +lammps_quest MD with quantum forces, coupling to Quest DFT code +lammps_spparks grain-growth Monte Carlo with strain via MD, + coupling to SPPARKS kinetic MC code +library collection of useful inter-code communication routines +simple simple example of driver code calling LAMMPS as library + +-------------------------------- + +The "simple" directory has its own README on how to proceed. The +driver code is provided both in C++ and in C. It simply runs a LAMMPS +simulation, extracts atom coordinates, changes a coordinate, passes +the coordinates back to LAMMPS, and runs some more. + +-------------------------------- + +The "library" directory has a small collection of routines, useful for +exchanging data between 2 codes being run together as a coupled +application. It is used by the LAMMPS <-> Quest and LAMMPS <-> +SPPARKS applications in the other 2 directories. + +The library dir has a Makefile (which you may need to edit for your +box). If you type + +g++ -f Makefile.g++ + +you should create libcouple.a, which the other coupled applications +link to. + +Note that the library uses MPI, so the Makefile you use needs to +include a path to the MPI include file, if it is not someplace +the compiler will find it. + +-------------------------------- + +The "lammps_quest" directory has an application that runs classical MD +via LAMMPS, but uses quantum forces calculated by the Quest DFT +(density functional) code in place of the usual classical MD forces +calculated by a pair style in LAMMPS. + +lmpqst.cpp main program + it links LAMMPS as a library + it invokes Quest as an executable +in.lammps LAMMPS input script, without the run command +si_111.in Quest input script for an 8-atom Si unit cell +lmppath.h contains path to LAMMPS home directory +qstexe.h contains full pathname to Quest executable + +After editing the Makefile, lmppath.h, and qstexe.h to make them +suitable for your box, type: + +g++ -f Makefile.g++ + +and you should get the lmpqst executable. + +NOTE: To run this coupled application, you must of course, have Quest +built on your system. It's WWW site is http://dft.sandia.gov/Quest. +It is not an open-source code, buy you can contact its authors to +obtain a copy. + +You can run lmpqst in serial or parallel as: + +% lmpqst Niter in.lammps in.quest +% mpirun -np 4 lmpqst Niter in.lammps in.quest + +where + +Niter = # of MD iterations +in.lammps = LAMMPS input script +in.quest = Quest input script + +The log files are for this run: + +% lmpqst 10 in.lammps si_111.in + +This application is an example of a coupling where the driver code +(lmpqst) runs one code (LAMMPS) as an outer code and facilitates it +calling the other code (Quest) as an inner code. Specifically, the +driver (lmpqst) invokes one code (LAMMPS) to perform its timestep +loop, and grabs information from the other code (Quest) during its +timestep. This is done in LAMMPS using the fix external command, +which makes a "callback" to the driver application (lmpqst), which in +turn invokes Quest with new atom coordinates, lets Quest compute +forces, and returns those forces to the LAMMPS fix external. + +The driver code launches LAMMPS in parallel. But Quest is only run on +a single processor. It would be possible to change this by using a +parallel build of Quest. + +Since Quest does not currently have a library interface, the driver +code interfaces with Quest via input and output files. + +Note that essentially 100% of the run time for this coupled +application is spent in Quest, as the quantum calculation of forces +dominates the calculation. + +You can look at the log files in the directory to see sample LAMMPS +output for this simulation. Dump files produced by LAMMPS are stored +as dump.md. + +-------------------------------- + +The "lammps_spparks" directory has an application that models grain +growth in the presence of strain. The grain growth is simulated by a +Potts model in a kinetic Monte Carlo code SPPARKS. Clusters of like +spins on a lattice represent grains. The Hamiltonian for the energy +due of a collection of spins includes a strain term and is described +on this page in the SPPARKS documentation: +http://www.sandia.gov/~sjplimp/spparks/doc/app_potts_strain.html. The +strain is computed by LAMMPS as a particle displacement where pairs of +atoms across a grain boundary are of different types and thus push off +from each other due to a Lennard-Jones sigma between particles of +different types that is larger than the sigma between particles of the +same type (interior to grains). + +lmpspk.cpp main program + it links LAMMPS and SPPARKS as libraries +in.spparks SPPARKS input script, without the run command +lmppath.h contains path to LAMMPS home directory +spkpath.h contains path to SPPARKS home directory + +After editing the Makefile, lmppath.h, and spkpath.h to make them +suitable for your box, type: + +g++ -f Makefile.g++ + +and you should get the lmpspk executable. + +NOTE: To build and run this coupled application, you must of course, +have SPPARKS built on your system. It's WWW site is +http://www.sandia.gov/~sjplimp/spparks.html. It is an open-source +code, written by two of the LAMMPS authors. + +You can run lmpspk in serial or parallel as: + +% lmpspk Niter Ndelta Sfactor in.spparks +% mpirun -np 4 lmpspk Niter Ndelta Sfactor in.spparks + +where + +Niter = # of outer iterations +Ndelta = time to run MC in each iteration +Sfactor = multiplier on strain effect +in.spparks = SPPARKS input script + +The log files are for this run: + +% lmpspk 5 10.0 5 in.spparks + +This application is an example of a coupling where the driver code +(lmpspk) alternates back and forth between the 2 applications (LAMMPS +and SPPARKS). Each outer timestep in the driver code, the following +tasks are performed. One code (SPPARKS) is invoked for a few Monte +Carlo steps. Some of its output (spin state) is passed to the other +code (LAMMPS) as input (atom type). The the other code (LAMMPS) is +invoked for a few timesteps. Some of its output (atom coords) is +massaged to become an input (per-atom strain) for the original code +(SPPARKS). + +The driver code launches both SPPARKS and LAMMPS in parallel and they +both decompose their spatial domains in the same manner. The datums +in SPPARKS (lattice sites) are the same as the datums in LAMMPS +(coarse-grained particles). If this were not the case, more +sophisticated inter-code communication could be performed. + +You can look at the log files in the directory to see sample LAMMPS +and SPPARKS output for this simulation. Dump files produced by the +run are stored as dump.mc and dump.md. The image*.png files show +snapshots from both the LAMMPS and SPPARKS output. Note that the +in.lammps and data.lammps files are not inputs; they are generated by +the lmpspk driver. diff --git a/couple/lammps_quest/Makefile.g++ b/couple/lammps_quest/Makefile.g++ new file mode 100644 index 0000000000..e9e4020091 --- /dev/null +++ b/couple/lammps_quest/Makefile.g++ @@ -0,0 +1,47 @@ +# Makefile for MD with quantum forces via LAMMPS <-> Quest coupling + +SHELL = /bin/sh + +# System-specific settings + +LAMMPS = /home/sjplimp/lammps + +CC = g++ +CCFLAGS = -g -O -DMPICH_IGNORE_CXX_SEEK -I../library +DEPFLAGS = -M +LINK = g++ +LINKFLAGS = -g -O -L../library -L${LAMMPS}/src +USRLIB = -lcouple -llmp_g++ +SYSLIB = -lfftw -lmpich -lpthread +ARCHIVE = ar +ARFLAGS = -rc +SIZE = size + +# Files + +EXE = lmpqst +SRC = $(wildcard *.cpp) +INC = $(wildcard *.h) +OBJ = $(SRC:.cpp=.o) + +# Targets + +$(EXE): $(OBJ) + $(LINK) $(LINKFLAGS) $(OBJ) $(USRLIB) $(SYSLIB) -o $(EXE) + $(SIZE) $(EXE) + +clean: + rm $(EXE) *.o + +# Compilation rules + +%.o:%.cpp + $(CC) $(CCFLAGS) -c $< + +%.d:%.cpp + $(CC) $(CCFLAGS) $(DEPFLAGS) $< > $@ + +# Individual dependencies + +DEPENDS = $(OBJ:.o=.d) +include $(DEPENDS) diff --git a/couple/lammps_quest/in.lammps b/couple/lammps_quest/in.lammps new file mode 100644 index 0000000000..c8d7faad6f --- /dev/null +++ b/couple/lammps_quest/in.lammps @@ -0,0 +1,20 @@ +# LAMMPS input for coupling MD/Quantum + +units metal +dimension 3 +atom_style atomic +atom_modify sort 0 0.0 + +lattice diamond 5.43 +region box block 0 1 0 1 0 1 +create_box 1 box +create_atoms 1 box +mass 1 28.08 + +velocity all create 300.0 87293 loop geom + +fix 1 all nve +fix 2 all external + +dump 1 all custom 1 dump.md id type x y z fx fy fz +thermo 1 diff --git a/couple/lammps_quest/lmppath.h b/couple/lammps_quest/lmppath.h new file mode 100644 index 0000000000..73f304abce --- /dev/null +++ b/couple/lammps_quest/lmppath.h @@ -0,0 +1 @@ +#define LMPPATH /home/sjplimp/lammps diff --git a/couple/lammps_quest/lmpqst.cpp b/couple/lammps_quest/lmpqst.cpp new file mode 100644 index 0000000000..661d8c04cd --- /dev/null +++ b/couple/lammps_quest/lmpqst.cpp @@ -0,0 +1,262 @@ +// lmpqst = umbrella driver to couple LAMMPS + Quest +// for MD using quantum forces + +// Syntax: lmpqst Niter in.lammps in.quest +// Niter = # of MD iterations +// in.lammps = LAMMPS input script +// in.quest = Quest input script + +#include "mpi.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" + +#include "many2one.h" +#include "one2many.h" +#include "files.h" +#include "memory.h" +#include "error.h" + +#define QUOTE_(x) #x +#define QUOTE(x) QUOTE_(x) + +#include "lmppath.h" +#include QUOTE(LMPPATH/src/lammps.h) +#include QUOTE(LMPPATH/src/library.h) +#include QUOTE(LMPPATH/src/input.h) +#include QUOTE(LMPPATH/src/modify.h) +#include QUOTE(LMPPATH/src/fix.h) +#include QUOTE(LMPPATH/src/fix_external.h) + +#include "qstexe.h" + +using namespace LAMMPS_NS; + +#define ANGSTROM_per_BOHR 0.529 +#define EV_per_RYDBERG 13.6056923 + +void quest_callback(void *, int, int, int *, double **, double **); + +struct Info { + int me; + Memory *memory; + LAMMPS *lmp; + char *quest_input; +}; + +/* ---------------------------------------------------------------------- */ + +int main(int narg, char **arg) +{ + int n; + char str[128]; + + // setup MPI + + MPI_Init(&narg,&arg); + MPI_Comm comm = MPI_COMM_WORLD; + + int me,nprocs; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + Memory *memory = new Memory(comm); + Error *error = new Error(comm); + + // command-line args + + if (narg != 4) error->all("Syntax: lmpqst Niter in.lammps in.quest"); + + int niter = atoi(arg[1]); + n = strlen(arg[2]) + 1; + char *lammps_input = new char[n]; + strcpy(lammps_input,arg[2]); + n = strlen(arg[3]) + 1; + char *quest_input = new char[n]; + strcpy(quest_input,arg[3]); + + // instantiate LAMMPS + + LAMMPS *lmp = new LAMMPS(0,NULL,MPI_COMM_WORLD); + + // create simulation in LAMMPS from in.lammps + + lmp->input->file(lammps_input); + + // make info avaiable to callback function + + Info info; + info.me = me; + info.memory = memory; + info.lmp = lmp; + info.quest_input = quest_input; + + // set callback to Quest inside fix external + + int ifix = lmp->modify->find_fix("2"); + FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix]; + fix->set_callback(quest_callback,&info); + + // run LAMMPS for Niter + // each time it needs forces, it will invoke quest_callback + + sprintf(str,"run %d",niter); + lmp->input->one(str); + + // clean up + + delete lmp; + + delete memory; + delete error; + + delete [] lammps_input; + delete [] quest_input; + + MPI_Finalize(); +} + +/* ---------------------------------------------------------------------- + callback to Quest with atom IDs and coords from each proc + invoke Quest to compute forces, load them into f for LAMMPS to use + f can be NULL if proc owns no atoms +------------------------------------------------------------------------- */ + +void quest_callback(void *ptr, int ntimestep, int nlocal, + int *id, double **x, double **f) +{ + int i,j; + char str[128]; + + Info *info = (Info *) ptr; + + // boxlines = LAMMPS box size converted into Quest lattice vectors + + char **boxlines = NULL; + if (info->me == 0) { + boxlines = new char*[3]; + for (i = 0; i < 3; i++) boxlines[i] = new char[128]; + } + + double boxxlo = *((double *) lammps_extract_global(info->lmp,"boxxlo")); + double boxxhi = *((double *) lammps_extract_global(info->lmp,"boxxhi")); + double boxylo = *((double *) lammps_extract_global(info->lmp,"boxylo")); + double boxyhi = *((double *) lammps_extract_global(info->lmp,"boxyhi")); + double boxzlo = *((double *) lammps_extract_global(info->lmp,"boxzlo")); + double boxzhi = *((double *) lammps_extract_global(info->lmp,"boxzhi")); + + double xprd = (boxxhi-boxxlo)/ANGSTROM_per_BOHR; + double yprd = (boxyhi-boxylo)/ANGSTROM_per_BOHR; + double zprd = (boxzhi-boxzlo)/ANGSTROM_per_BOHR; + + if (info->me == 0) { + sprintf(boxlines[0],"%g %g %g\n",xprd,0.0,0.0); + sprintf(boxlines[1],"%g %g %g\n",0.0,yprd,0.0); + sprintf(boxlines[2],"%g %g %g\n",0.0,0.0,zprd); + } + + // xlines = x for atoms on each proc converted to text lines + // xlines is suitable for insertion into Quest input file + // convert LAMMPS Angstroms to Quest bohr + + int natoms; + MPI_Allreduce(&nlocal,&natoms,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); + + Many2One *lmp2qst = new Many2One(MPI_COMM_WORLD); + lmp2qst->setup(nlocal,id,natoms); + + char **xlines = NULL; + double **xquest = NULL; + if (info->me == 0) { + xquest = info->memory->create_2d_double_array(natoms,3,"lmpqst:xquest"); + xlines = new char*[natoms]; + for (i = 0; i < natoms; i++) xlines[i] = new char[128]; + } + + if (info->me == 0) lmp2qst->gather(&x[0][0],3,&xquest[0][0]); + else lmp2qst->gather(&x[0][0],3,NULL); + + if (info->me == 0) { + for (i = 0; i < natoms; i++) { + xquest[i][0] /= ANGSTROM_per_BOHR; + xquest[i][1] /= ANGSTROM_per_BOHR; + xquest[i][2] /= ANGSTROM_per_BOHR; + } + for (i = 0; i < natoms; i++) { + sprintf(xlines[i],"%d %d %g %g %g\n",i+1,1, + xquest[i][0],xquest[i][1],xquest[i][2]); + } + } + + // one-processor tasks: + // whack all lcao.* files + // cp quest_input to lcao.in + // replace atom coords section of lcao.in with new atom coords + // run Quest on one proc, save screen output to file + // flines = atom forces extracted from Quest screen file + // fquest = atom forces + // convert Quest Ryd/bohr to LAMMPS eV/Angstrom + + char **flines = NULL; + double **fquest = NULL; + if (info->me == 0) { + fquest = info->memory->create_2d_double_array(natoms,3,"lmpqst:fquest"); + flines = new char*[natoms]; + for (i = 0; i < natoms; i++) flines[i] = new char[128]; + } + + if (info->me == 0) { + system("rm lcao.*"); + sprintf(str,"cp %s lcao.in",info->quest_input); + system(str); + sprintf(str,"cp %s lcao.x",QUOTE(QUEST)); + system(str); + replace("lcao.in","primitive lattice vectors",3,boxlines); + replace("lcao.in","atom, type, position vector",natoms,xlines); + system("lcao.x > lcao.screen"); + extract("lcao.screen","atom x force " + "y force z force",natoms,flines); + + int itmp; + for (i = 0; i < natoms; i++) + sscanf(flines[i],"%d %lg %lg %lg",&itmp, + &fquest[i][0],&fquest[i][1],&fquest[i][2]); + + for (i = 0; i < natoms; i++) { + fquest[i][0] *= EV_per_RYDBERG / ANGSTROM_per_BOHR; + fquest[i][1] *= EV_per_RYDBERG / ANGSTROM_per_BOHR; + fquest[i][2] *= EV_per_RYDBERG / ANGSTROM_per_BOHR; + } + } + + // convert fquest on one proc into f for atoms on each proc + + One2Many *qst2lmp = new One2Many(MPI_COMM_WORLD); + qst2lmp->setup(natoms,nlocal,id); + double *fvec = NULL; + if (f) fvec = &f[0][0]; + if (info->me == 0) qst2lmp->scatter(&fquest[0][0],3,fvec); + else qst2lmp->scatter(NULL,3,fvec); + + // clean up + // some data only exists on proc 0 + + delete lmp2qst; + delete qst2lmp; + + info->memory->destroy_2d_double_array(xquest); + info->memory->destroy_2d_double_array(fquest); + + if (boxlines) { + for (i = 0; i < 3; i++) delete [] boxlines[i]; + delete [] boxlines; + } + if (xlines) { + for (i = 0; i < natoms; i++) delete [] xlines[i]; + delete [] xlines; + } + if (flines) { + for (i = 0; i < natoms; i++) delete [] flines[i]; + delete [] flines; + } +} diff --git a/couple/lammps_quest/qstexe.h b/couple/lammps_quest/qstexe.h new file mode 100644 index 0000000000..19d49bd7d9 --- /dev/null +++ b/couple/lammps_quest/qstexe.h @@ -0,0 +1 @@ +#define QUEST /home/sjplimp/csrf/quest/src/lcao.x diff --git a/couple/lammps_quest/si_111.in b/couple/lammps_quest/si_111.in new file mode 100644 index 0000000000..7c7eedfa48 --- /dev/null +++ b/couple/lammps_quest/si_111.in @@ -0,0 +1,161 @@ +do setup +do iters +do force +no relax +setup data +title +Si 1x1x1 unit cell +functional +PBE +dimensions of system (0=cluster ... 3=bulk) +3 +primitive lattice vectors +10.261212 0.000000 0.000000 +0.000000 10.261212 0.000000 +0.000000 0.000000 10.261212 +grid dimensions +10 10 10 +atom types +1 +type number, label: + 1 Si_pbe +notes5 + Originally constructed by Peter A. Schultz, 12Apr01 + potential generated by new Hamann program PUNSLDX + Cite use with: D.R. Hamann, unpublished. + Potential: "standard" setting out to l=2 + Basis: amended Jun05 for better (2d/1d not 1d/1d) d-function +effective nuclear charge (s2p2 to 10.0) + 4.00000000d+00 +pseudopotentials: Lmax, and effective gaussian range + 2 0.86000000d+00 +functional type used in generating potential: +PBE +radial mesh: number of points for local and non-local pot integrals + 80 67 +mesh points for nuclear potential; ham2dh + 0.02500000 0.02696978 0.02909477 0.03138719 0.03386023 0.03652812 + 0.03940622 0.04251109 0.04586060 0.04947402 0.05337215 0.05757741 + 0.06211402 0.06700807 0.07228773 0.07798338 0.08412779 0.09075634 + 0.09790716 0.10562140 0.11394345 0.12292121 0.13260635 0.14305458 + 0.15432605 0.16648562 0.17960325 0.19375443 0.20902061 0.22548964 + 0.24325628 0.26242278 0.28309943 0.30540522 0.32946852 0.35542780 + 0.38343245 0.41364362 0.44623518 0.48139466 0.51932441 0.56024270 + 0.60438500 0.65200533 0.70337773 0.75879783 0.81858456 0.88308197 + 0.95266121 1.02772271 1.10869840 1.19605428 1.29029305 1.39195702 + 1.50163124 1.61994684 1.74758469 1.88527930 2.03382306 2.19407079 + 2.36694466 2.55343950 2.75462852 2.97166951 3.20581145 3.45840177 + 3.73089402 4.02485632 4.34198031 4.68409093 5.05315693 5.45130215 + 5.88081777 6.34417553 6.84404189 7.38329340 7.96503329 8.59260927 + 9.26963282 10.00000000 +radwts: weights for radial points + 0.00189603 0.00204542 0.00220659 0.00238045 0.00256800 0.00277034 + 0.00298862 0.00322410 0.00347813 0.00375218 0.00404781 0.00436675 + 0.00471081 0.00508198 0.00548240 0.00591436 0.00638036 0.00688308 + 0.00742541 0.00801047 0.00864162 0.00932251 0.01005704 0.01084945 + 0.01170429 0.01262649 0.01362135 0.01469459 0.01585240 0.01710143 + 0.01844888 0.01990249 0.02147064 0.02316234 0.02498733 0.02695611 + 0.02908002 0.03137128 0.03384307 0.03650961 0.03938625 0.04248955 + 0.04583736 0.04944895 0.05334510 0.05754823 0.06208254 0.06697411 + 0.07225109 0.07794385 0.08408515 0.09071034 0.09785753 0.10556786 + 0.11388570 0.12285891 0.13253914 0.14298208 0.15424783 0.16640123 + 0.17951222 0.19365623 0.20891467 0.22537535 0.24313298 0.26228977 + 0.28295594 0.30525043 0.32930153 0.35524766 0.38323811 0.41343397 + 0.44600900 0.48115067 0.51906119 0.55995874 0.60407867 0.65167486 + 0.70302122 0.75841323 +non-local potential: l,potential*integration weight + 0 0.62022930 0.62128855 0.62243016 0.62366033 0.62498568 0.62641328 + 0.62795061 0.62960563 0.63138673 0.63330275 0.63536294 0.63757692 + 0.63995464 0.64250630 0.64524218 0.64817253 0.65130735 0.65465605 + 0.65822713 0.66202767 0.66606269 0.67033437 0.67484108 0.67957602 + 0.68452576 0.68966817 0.69497006 0.70038419 0.70584566 0.71126756 + 0.71653578 0.72150290 0.72598113 0.72973436 0.73246932 0.73382636 + 0.73337030 0.73058243 0.72485505 0.71549107 0.70171167 0.68267654 + 0.65752236 0.62542611 0.58570073 0.53792896 0.48213811 0.41900888 + 0.35009536 0.27800640 0.20646172 0.14009458 0.08384960 0.04186877 + 0.01596164 0.00423035 0.00115036 0.00066636 0.00047879 0.00029939 + 0.00016329 0.00007995 0.00003517 0.00001362 0.00000445 0.00000111 + 0.00000016 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 +non-local potential: l,potential*integration weight + 1 0.59551624 0.59463303 0.59368033 0.59265268 0.59154422 0.59034862 + 0.58905906 0.58766819 0.58616811 0.58455033 0.58280567 0.58092430 + 0.57889565 0.57670833 0.57435015 0.57180802 0.56906791 0.56611482 + 0.56293268 0.55950435 0.55581158 0.55183493 0.54755377 0.54294628 + 0.53798942 0.53265896 0.52692951 0.52077458 0.51416671 0.50707751 + 0.49947790 0.49133817 0.48262822 0.47331766 0.46337588 0.45277197 + 0.44147437 0.42945016 0.41666374 0.40307468 0.38863443 0.37328165 + 0.35693601 0.33949042 0.32080256 0.30068740 0.27891443 0.25521609 + 0.22931791 0.20100526 0.17024474 0.13737521 0.10336405 0.07007167 + 0.04035673 0.01767907 0.00470635 0.00076638 0.00047880 0.00029939 + 0.00016329 0.00007995 0.00003517 0.00001362 0.00000445 0.00000111 + 0.00000016 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 +non-local potential: l,potential*integration weight + 2 0.56305372 0.55961728 0.55591134 0.55191498 0.54760572 0.54295941 + 0.53795013 0.53255008 0.52672947 0.52045641 0.51369682 0.50641433 + 0.49857022 0.49012333 0.48103004 0.47124429 0.46071759 0.44939919 + 0.43723624 0.42417413 0.41015690 0.39512792 0.37903070 0.36181001 + 0.34341340 0.32379300 0.30290805 0.28072780 0.25723539 0.23243242 + 0.20634465 0.17902876 0.15058041 0.12114359 0.09092117 0.06018665 + 0.02929636 -0.00129833 -0.03104046 -0.05926034 -0.08517498 -0.10789810 + -0.12646610 -0.13988656 -0.14721657 -0.14767751 -0.14080976 -0.12666296 + -0.10600305 -0.08049270 -0.05276798 -0.02629475 -0.00486427 0.00837657 + 0.01228139 0.00892332 0.00342796 0.00074936 0.00047880 0.00029939 + 0.00016329 0.00007995 0.00003517 0.00001362 0.00000445 0.00000111 + 0.00000016 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 +number of radial functions **** Si PBE Ham-II basis 20Feb01-PAS **** + 5 +angular momentum, number of alphas + 0 4 +alphas - s - 4s/2/4s407 (bulk Si dzp Eopt+reopt c1) + 0.10460000d+00 0.27226300d+00 1.30050800d+00 2.60103000d+00 +wave function coefficients + 0.20995300d+00 0.55978200d+00 -0.99128200d+00 0.33487100d+00 +angular momentum, number of alphas + 1 3 +alphas - p - 3p/2/3p492 (bulk Si dzp Eopt + reopt c1) + 0.09424100d+00 0.31767900d+00 1.56114500d+00 +wave function coefficients + 0.06761600d+00 0.31821200d+00 -0.06638300d+00 +angular momentum, number of alphas + 0 1 +alphas - s - second zeta s polarization + 0.10460000d+00 +wave function coefficients + 1.00000000d+00 +angular momentum, number of alphas + 1 1 +alphas - p - second zeta p polarization + 0.09424100d+00 +wave function coefficients + 1.00000000d+00 +angular momentum, number of alphas + 2 2 +alphas - d - angular polarization (dzp Eopt) + 0.32000000d+00 1.40000000d+00 +wave function coefficients + 0.31557000d+00 1.00000000d+00 +shell occupancies for this silicon, Si: s(2.00)p(2.00) + 2.00000000 2.00000000 0.00000000 0.00000000 0.00000000 0.00000000 +end atom file +number of atoms in unit cell +8 +atom, type, position vector +1 1 0.0000000000 0.0000000000 0.0000000000 +2 1 5.1306060590 5.1306060590 0.0000000000 +3 1 5.1306060590 0.0000000000 5.1306060590 +4 1 0.0000000000 5.1306060590 5.1306060590 +5 1 2.5653030295 2.5653030295 2.5653030295 +6 1 7.6959090885 7.6959090885 2.5653030295 +7 1 7.6959090885 2.5653030295 7.6959090885 +8 1 2.5653030295 7.6959090885 7.6959090885 +kgrid +0 0 0 +end setup phase data +run phase input data +end of run phase data diff --git a/couple/lammps_spparks/Makefile b/couple/lammps_spparks/Makefile new file mode 100644 index 0000000000..df90a14e4a --- /dev/null +++ b/couple/lammps_spparks/Makefile @@ -0,0 +1,48 @@ +# Makefile for grain growth via LAMMPS <-> SPPARKS coupling + +SHELL = /bin/sh + +# System-specific settings + +LAMMPS = /home/sjplimp/lammps +SPPARKS = /home/sjplimp/spparks + +CC = g++ +CCFLAGS = -g -O -DMPICH_IGNORE_CXX_SEEK -I../library +DEPFLAGS = -M +LINK = g++ +LINKFLAGS = -g -O -L../library -L${LAMMPS}/src -L${SPPARKS}/src +USRLIB = -lcouple -llmp_g++ -lspk_g++ +SYSLIB = -lfftw -lmpich -lpthread +ARCHIVE = ar +ARFLAGS = -rc +SIZE = size + +# Files + +EXE = lmpspk +SRC = $(wildcard *.cpp) +INC = $(wildcard *.h) +OBJ = $(SRC:.cpp=.o) + +# Targets + +lmpspk: $(OBJ) + $(LINK) $(LINKFLAGS) $(OBJ) $(USRLIB) $(SYSLIB) -o $(EXE) + $(SIZE) $(EXE) + +clean: + rm $(EXE) *.o + +# Compilation rules + +%.o:%.cpp + $(CC) $(CCFLAGS) -c $< + +%.d:%.cpp + $(CC) $(CCFLAGS) $(DEPFLAGS) $< > $@ + +# Individual dependencies + +DEPENDS = $(OBJ:.o=.d) +include $(DEPENDS) diff --git a/couple/lammps_spparks/in.spparks b/couple/lammps_spparks/in.spparks new file mode 100644 index 0000000000..c2bb5b74b2 --- /dev/null +++ b/couple/lammps_spparks/in.spparks @@ -0,0 +1,23 @@ +# SPPARKS input for coupling MD/MC + +seed 56789 + +app_style potts/strain 100 + +dimension 2 +lattice sq/8n 1.0 +region box block 0 50 0 50 -0.5 0.5 +create_box box +create_sites box +set site range 1 100 +set d1 value 0.0 + +sector yes +solve_style tree + +diag_style energy + +temperature 1.0 + +stats 10.0 +dump 1 10.0 dump.mc diff --git a/couple/lammps_spparks/lmppath.h b/couple/lammps_spparks/lmppath.h new file mode 100644 index 0000000000..73f304abce --- /dev/null +++ b/couple/lammps_spparks/lmppath.h @@ -0,0 +1 @@ +#define LMPPATH /home/sjplimp/lammps diff --git a/couple/lammps_spparks/lmpspk.cpp b/couple/lammps_spparks/lmpspk.cpp new file mode 100644 index 0000000000..da6d252369 --- /dev/null +++ b/couple/lammps_spparks/lmpspk.cpp @@ -0,0 +1,223 @@ +// lmpspk = umbrella driver to couple LAMMPS + SPPARKS +// for a strain-induced grain growth model + +// Syntax: lmpspk Niter Ndelta Sfactor in.spparks +// Niter = # of outer iterations +// Ndelta = time to run MC in each iteration +// Sfactor = multiplier on strain effect +// in.spparks = SPPARKS input script + +#include "mpi.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" + +#include "lammps_data_write.h" +#include "many2many.h" +#include "memory.h" +#include "error.h" + +#define QUOTE_(x) #x +#define QUOTE(x) QUOTE_(x) + +#include "spkpath.h" +#include QUOTE(SPKPATH/src/spparks.h) +#include QUOTE(SPKPATH/src/library.h) +#include QUOTE(SPKPATH/src/input.h) + +#include "lmppath.h" +#include QUOTE(LMPPATH/src/lammps.h) +#include QUOTE(LMPPATH/src/library.h) +#include QUOTE(LMPPATH/src/input.h) +#include QUOTE(LMPPATH/src/modify.h) +#include QUOTE(LMPPATH/src/compute.h) + +using namespace SPPARKS_NS; +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +int main(int narg, char **arg) +{ + int i,n; + char str[128]; + + // setup MPI + + MPI_Init(&narg,&arg); + MPI_Comm comm = MPI_COMM_WORLD; + + int me,nprocs; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + Memory *memory = new Memory(comm); + Error *error = new Error(comm); + + // command-line args + + if (narg != 5) + error->all("Syntax: lmpspk Niter Ndelta Sfactor in.spparks"); + + int niter = atoi(arg[1]); + double delta = atof(arg[2]); + double sfactor = atof(arg[3]); + n = strlen(arg[4]) + 1; + char *spparks_input = new char[n]; + strcpy(spparks_input,arg[4]); + + // instantiate LAMMPS and SPPARKS + + SPPARKS *spk = new SPPARKS(0,NULL,MPI_COMM_WORLD); + LAMMPS *lmp = new LAMMPS(0,NULL,MPI_COMM_WORLD); + + // create simulation in SPPARKS from in.spparks + + spk->input->file(spparks_input); + + // extract permanent info from SPPARKS + + int dimension,nglobal,nlocal_spparks,nspins; + double boxxlo,boxxhi,boxylo,boxyhi,boxzlo,boxzhi; + int *id_spparks,*spins; + double **xyz; + double *strain; + + dimension = *((int *) spparks_extract(spk,"dimension")); + nglobal = *((int *) spparks_extract(spk,"nglobal")); + nlocal_spparks = *((int *) spparks_extract(spk,"nlocal")); + + boxxlo = *((double *) spparks_extract(spk,"boxxlo")); + boxxhi = *((double *) spparks_extract(spk,"boxxhi")); + boxylo = *((double *) spparks_extract(spk,"boxylo")); + boxyhi = *((double *) spparks_extract(spk,"boxyhi")); + if (dimension == 3) { + boxzlo = *((double *) spparks_extract(spk,"boxzlo")); + boxzhi = *((double *) spparks_extract(spk,"boxzhi")); + } else { + boxzlo = -0.5; + boxzhi = 0.5; + } + + id_spparks = (int *) spparks_extract(spk,"id"); + spins = (int *) spparks_extract(spk,"site"); + xyz = (double **) spparks_extract(spk,"xyz"); + + nspins = *((int *) spparks_extract(spk,"nspins")); + strain = (double *) spparks_extract(spk,"strain"); + + // write a LAMMPS input script using SPPARKS params + + if (me == 0) { + FILE *fp = fopen("in.lammps","w"); + if (fp == NULL) error->one("Could not create LAMMPS input script"); + + fprintf(fp,"units lj\n"); + sprintf(str,"dimension %d\n",dimension); + fprintf(fp,str); + fprintf(fp,"atom_style atomic\n\n"); + + fprintf(fp,"read_data data.lammps\n"); + fprintf(fp,"mass * 1.0\n\n"); + + fprintf(fp,"pair_style lj/cut 2.5\n"); + fprintf(fp,"pair_coeff * * 1.0 1.2\n"); + for (i = 0; i < nspins; i++) { + sprintf(str,"pair_coeff %d %d 1.0 1.0\n",i+1,i+1); + fprintf(fp,str); + } + fprintf(fp,"\n"); + + fprintf(fp,"compute da all displace/atom\n\n"); + fprintf(fp,"dump 1 all atom 10 dump.md\n"); + fprintf(fp,"thermo 1\n"); + + fclose(fp); + } + + // write a LAMMPS data file using SPPARKS data + + LAMMPSDataWrite *lwd = new LAMMPSDataWrite(MPI_COMM_WORLD); + lwd->file("data.lammps"); + lwd->header("%d atoms",nglobal); + lwd->header("%d atom types",nspins); + lwd->header("%g %g xlo xhi",boxxlo,boxxhi); + lwd->header("%g %g ylo yhi",boxylo,boxyhi); + lwd->header("%g %g zlo zhi",boxzlo,boxzhi); + lwd->atoms(nlocal_spparks); + lwd->atoms(id_spparks); + lwd->atoms(spins); + lwd->atoms(3,xyz); + lwd->execute(); + delete lwd; + + // create simulation in LAMMPS from created input script + + lmp->input->file("in.lammps"); + + // create transfer operators + + Many2Many *spk2lmp = new Many2Many(MPI_COMM_WORLD); + Many2Many *lmp2spk = new Many2Many(MPI_COMM_WORLD); + + // timestep loop + // run SPPARKS for delta time + // use SPPARKS spins to reset LAMMPS atom types + // perform LAMMPS minimization + // use atom displacements to reset strain values in SPPARKS + // realloc displace as necessary since nlocal_lammps may change + // re-create both xfers every iteration since LAMMPS may migrate atoms + + int nmax = 0; + double *displace = NULL; + + int nlocal_lammps; + int *id_lammps,*type; + double **displace_lammps; + + for (int iter = 0; iter < niter; iter++) { + sprintf(str,"run %g",delta); + spk->input->one(str); + + nlocal_lammps = *((int *) lammps_extract_global(lmp,"nlocal")); + id_lammps = (int *) lammps_extract_atom(lmp,"id"); + type = (int *) lammps_extract_atom(lmp,"type"); + + spk2lmp->setup(nlocal_spparks,id_spparks,nlocal_lammps,id_lammps); + spk2lmp->exchange(spins,type); + + lmp->input->one("minimize 0.001 0.001 10 1000"); + + nlocal_lammps = *((int *) lammps_extract_global(lmp,"nlocal")); + id_lammps = (int *) lammps_extract_atom(lmp,"id"); + displace_lammps = (double **) lammps_extract_compute(lmp,"da",1,2); + + if (nlocal_lammps > nmax) { + memory->sfree(displace); + nmax = nlocal_lammps; + displace = (double *) + memory->smalloc(nmax*sizeof(double),"lmpspk:displace"); + } + + for (i = 0; i < nlocal_lammps; i++) + displace[i] = sfactor*displace_lammps[i][3]; + + lmp2spk->setup(nlocal_lammps,id_lammps,nlocal_spparks,id_spparks); + lmp2spk->exchange(displace,strain); + } + + memory->sfree(displace); + + // clean up + + delete spk2lmp; + delete lmp2spk; + delete spk; + delete lmp; + + delete [] spparks_input; + delete memory; + delete error; + + MPI_Finalize(); +} diff --git a/couple/lammps_spparks/spkpath.h b/couple/lammps_spparks/spkpath.h new file mode 100644 index 0000000000..2c1d959829 --- /dev/null +++ b/couple/lammps_spparks/spkpath.h @@ -0,0 +1 @@ +#define SPKPATH /home/sjplimp/spparks diff --git a/couple/library/Makefile.g++ b/couple/library/Makefile.g++ new file mode 100644 index 0000000000..07e5807680 --- /dev/null +++ b/couple/library/Makefile.g++ @@ -0,0 +1,42 @@ +# Makefile for coupling library + +SHELL = /bin/sh + +# System-specific settings + +CC = g++ +CCFLAGS = -g -O -DMPICH_IGNORE_CXX_SEEK +DEPFLAGS = -M +LINK = g++ +LINKFLAGS = -g -O +ARCHIVE = ar +ARFLAGS = -rc +SIZE = size + +# Files + +LIB = libcouple.a +SRC = $(wildcard *.cpp) +INC = $(wildcard *.h) +OBJ = $(SRC:.cpp=.o) + +# Targets + +lib: $(OBJ) + $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + +clean: + rm $(LIB) *.o *.d + +# Compilation rules + +%.o:%.cpp + $(CC) $(CCFLAGS) -c $< + +%.d:%.cpp + $(CC) $(CCFLAGS) $(DEPFLAGS) $< > $@ + +# Individual dependencies + +DEPENDS = $(OBJ:.o=.d) +include $(DEPENDS) diff --git a/couple/library/error.cpp b/couple/library/error.cpp new file mode 100644 index 0000000000..c70ec26f93 --- /dev/null +++ b/couple/library/error.cpp @@ -0,0 +1,42 @@ +#include "mpi.h" +#include "stdlib.h" +#include "stdio.h" +#include "error.h" + +/* ---------------------------------------------------------------------- */ + +Error::Error(MPI_Comm caller) +{ + comm = caller; + MPI_Comm_rank(comm,&me); +} + +/* ---------------------------------------------------------------------- + called by all procs +------------------------------------------------------------------------- */ + +void Error::all(const char *str) +{ + if (me == 0) printf("ERROR: %s\n",str); + MPI_Finalize(); + exit(1); +} + +/* ---------------------------------------------------------------------- + called by one proc +------------------------------------------------------------------------- */ + +void Error::one(const char *str) +{ + printf("ERROR on proc %d: %s\n",me,str); + MPI_Abort(comm,1); +} + +/* ---------------------------------------------------------------------- + called by one proc +------------------------------------------------------------------------- */ + +void Error::warning(const char *str) +{ + printf("WARNING: %s\n",str); +} diff --git a/couple/library/error.h b/couple/library/error.h new file mode 100644 index 0000000000..f62f620585 --- /dev/null +++ b/couple/library/error.h @@ -0,0 +1,19 @@ +#ifndef ERROR_H +#define ERROR_H + +#include "mpi.h" + +class Error { + public: + Error(MPI_Comm); + + void all(const char *); + void one(const char *); + void warning(const char *); + + private: + MPI_Comm comm; + int me; +}; + +#endif diff --git a/couple/library/files.cpp b/couple/library/files.cpp new file mode 100644 index 0000000000..83464deb32 --- /dev/null +++ b/couple/library/files.cpp @@ -0,0 +1,48 @@ +#include "stdio.h" +#include "string.h" +#include "files.h" + +#define MAXLINE 256 + +/* ---------------------------------------------------------------------- */ + +void replace(char *file, char *header, int n, char **lines) +{ + FILE *fpr = fopen(file,"r"); + FILE *fpw = fopen("tmp.file","w"); + + char line[MAXLINE]; + while (fgets(line,MAXLINE,fpr)) { + if (strstr(line,header)) { + fprintf(fpw,"%s",line); + for (int i = 0; i < n; i++) { + fgets(line,MAXLINE,fpr); + fprintf(fpw,"%s",lines[i]); + } + } else fprintf(fpw,"%s",line); + } + + fclose(fpr); + fclose(fpw); + rename("tmp.file",file); +} + +/* ---------------------------------------------------------------------- */ + +char **extract(char *file, char *header, int n, char **lines) +{ + FILE *fp = fopen(file,"r"); + + char line[MAXLINE]; + while (fgets(line,MAXLINE,fp)) { + if (strstr(line,header)) { + for (int i = 0; i < n; i++) { + fgets(line,MAXLINE,fp); + sprintf(lines[i],"%s",line); + } + break; + } + } + + fclose(fp); +} diff --git a/couple/library/files.h b/couple/library/files.h new file mode 100644 index 0000000000..75812d2991 --- /dev/null +++ b/couple/library/files.h @@ -0,0 +1,2 @@ +void replace(char *, char *, int, char **); +char **extract(char *, char *, int, char **); diff --git a/couple/library/irregular.cpp b/couple/library/irregular.cpp new file mode 100644 index 0000000000..aea9b8332c --- /dev/null +++ b/couple/library/irregular.cpp @@ -0,0 +1,393 @@ +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "irregular.h" +#include "memory.h" +#include "error.h" + +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +enum{UNSET,SET}; +enum{NONE,SAME,VARYING}; + +/* ---------------------------------------------------------------------- */ + +Irregular::Irregular(MPI_Comm caller) +{ + comm = caller; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + memory = new Memory(comm); + error = new Error(comm); + + init(); + + patternflag = UNSET; + sizestyle = NONE; +} + +/* ---------------------------------------------------------------------- */ + +Irregular::~Irregular() +{ + delete memory; + delete error; + deallocate(); +} + +/* ---------------------------------------------------------------------- + n = # of datums contributed by this proc + proclist = which proc each datum is to be sent to +------------------------------------------------------------------------- */ + +void Irregular::pattern(int n, int *proclist) +{ + // free any previous irregular info + + deallocate(); + init(); + + patternflag = SET; + sizestyle = NONE; + + ndatumsend = n; + + // list = 1 for procs I send to, including self + // nrecv = # of messages I receive, not including self + // self = 0 if no data for self, 1 if there is + + int *list = new int[nprocs]; + int *counts = new int[nprocs]; + + for (int i = 0; i < nprocs; i++) { + list[i] = 0; + counts[i] = 1; + } + + for (int i = 0; i < n; i++) list[proclist[i]] = 1; + MPI_Reduce_scatter(list,&nrecv,counts,MPI_INT,MPI_SUM,comm); + + self = 0; + if (list[me]) self = 1; + if (self) nrecv--; + + // storage for recv info, not including self + + recvproc = new int[nrecv]; + recvcount = new int[nrecv]; + recvsize = new int[nrecv]; + request = new MPI_Request[nrecv]; + status = new MPI_Status[nrecv]; + + // list = # of datums to send to each proc, including self + // nsend = # of messages I send, not including self + + for (int i = 0; i < nprocs; i++) list[i] = 0; + for (int i = 0; i < n; i++) list[proclist[i]]++; + + nsend = 0; + for (int i = 0; i < nprocs; i++) if (list[i] > 0) nsend++; + if (self) nsend--; + + // storage for send info, including self + + sendproc = new int[nsend+self]; + sendcount = new int[nsend+self]; + sendsize = new int[nsend+self]; + sendindices = (int *) memory->smalloc(n*sizeof(int),"sendindices"); + + // setup sendprocs and sendcounts, including self + // each proc begins with iproc > me, and continues until iproc = me + // list ends up with pointer to which send that proc is associated with + + int iproc = me; + int isend = 0; + for (int i = 0; i < nprocs; i++) { + iproc++; + if (iproc == nprocs) iproc = 0; + if (list[iproc] > 0) { + sendproc[isend] = iproc; + sendcount[isend] = list[iproc]; + list[iproc] = isend; + isend++; + } + } + + // post all receives for datum counts + + for (int i = 0; i < nrecv; i++) + MPI_Irecv(&recvcount[i],1,MPI_INT,MPI_ANY_SOURCE,0,comm,&request[i]); + + // barrier to insure receives are posted + + MPI_Barrier(comm); + + // send each datum count, packing buf with needed datums + + for (int i = 0; i < nsend; i++) + MPI_Send(&sendcount[i],1,MPI_INT,sendproc[i],0,comm); + + // insure all MPI_ANY_SOURCE messages are received + // set recvproc + + if (nrecv) MPI_Waitall(nrecv,request,status); + for (int i = 0; i < nrecv; i++) recvproc[i] = status[i].MPI_SOURCE; + + // ndatumrecv = total datums received, including self + + ndatumrecv = 0; + for (int i = 0; i < nrecv; i++) + ndatumrecv += recvcount[i]; + if (self) ndatumrecv += sendcount[nsend]; + + // setup sendindices, including self + // counts = offset into sendindices for each proc I send to + // let sc0 = sendcount[0], sc1 = sendcount[1], etc + // sendindices[0:sc0-1] = indices of datums in 1st message + // sendindices[sc0:sc0+sc1-1] = indices of datums in 2nd message, etc + + counts[0] = 0; + for (int i = 1; i < nsend+self; i++) + counts[i] = counts[i-1] + sendcount[i-1]; + + for (int i = 0; i < n; i++) { + isend = list[proclist[i]]; + sendindices[counts[isend]++] = i; + } + + // clean up + + delete [] counts; + delete [] list; +} + +/* ---------------------------------------------------------------------- + n = size of each received datum + return total size in bytes of received data on this proc +------------------------------------------------------------------------- */ + +int Irregular::size(int n) +{ + if (patternflag == UNSET) error->all("Cannot size without pattern"); + sizestyle = SAME; + + nsize = n; + + nsendmax = 0; + for (int i = 0; i < nsend+self; i++) { + sendsize[i] = nsize * sendcount[i]; + if (i < nsend) nsendmax = MAX(nsendmax,sendsize[i]); + } + + for (int i = 0; i < nrecv; i++) recvsize[i] = nsize * recvcount[i]; + nbytesrecv = nsize * ndatumrecv; + + return nbytesrecv; +} + +/* ---------------------------------------------------------------------- + slength,rlength = size of each datum to send and recv + soffset = offset into eventual buffer of send data for each datum + soffset can be NULL, in which case will build sendoffset from slength + return total size in bytes of received data on this proc +------------------------------------------------------------------------- */ + +int Irregular::size(int *slength, int *soffset, int *rlength) +{ + if (patternflag == UNSET) error->all("Cannot size without pattern"); + sizestyle = VARYING; + + // store local copy of pointers to send lengths/offsets + // if soffset not provided, create local copy from slength + + sendsizedatum = slength; + + if (soffset == NULL) { + sendoffsetflag = 1; + sendoffset = (int *) memory->smalloc(ndatumsend*sizeof(int),"sendoffset"); + + if (ndatumsend) sendoffset[0] = 0; + for (int i = 1; i < ndatumsend; i++) + sendoffset[i] = sendoffset[i-1] + sendsizedatum[i-1]; + + } else { + if (sendoffsetflag) memory->sfree(sendoffset); + sendoffsetflag = 0; + sendoffset = soffset; + } + + nsendmax = 0; + int m = 0; + for (int i = 0; i < nsend+self; i++) { + sendsize[i] = 0; + for (int j = 0; j < sendcount[i]; j++) + sendsize[i] += sendsizedatum[sendindices[m++]]; + if (i < nsend) nsendmax = MAX(nsendmax,sendsize[i]); + } + + nbytesrecv = 0; + m = 0; + for (int i = 0; i < nrecv; i++) { + recvsize[i] = 0; + for (int j = 0; j < recvcount[i]; j++) recvsize[i] += rlength[m++]; + nbytesrecv += recvsize[i]; + } + if (self) nbytesrecv += sendsize[nsend]; + + return nbytesrecv; +} + +/* ---------------------------------------------------------------------- + wrapper on 2 versions of exchange +------------------------------------------------------------------------- */ + +void Irregular::exchange(char *sendbuf, char *recvbuf) +{ + if (sizestyle == SAME) exchange_same(sendbuf,recvbuf); + else if (sizestyle == VARYING) exchange_varying(sendbuf,recvbuf); + else error->all("Irregular size was not set"); +} + +/* ---------------------------------------------------------------------- + sendbuf = data to send + recvbuf = buffer to recv all data into + requires nsize,nsendmax,recvsize,sendsize be setup by size(int) +------------------------------------------------------------------------- */ + +void Irregular::exchange_same(char *sendbuf, char *recvbuf) +{ + // post all receives + + int recvoffset = 0; + for (int irecv = 0; irecv < nrecv; irecv++) { + MPI_Irecv(&recvbuf[recvoffset],recvsize[irecv],MPI_BYTE, + recvproc[irecv],0,comm,&request[irecv]); + recvoffset += recvsize[irecv]; + } + + // malloc buf for largest send + + char *buf = (char *) memory->smalloc(nsendmax,"buf"); + + // barrier to insure receives are posted + + MPI_Barrier(comm); + + // send each message, packing buf with needed datums + + int m = 0; + for (int isend = 0; isend < nsend; isend++) { + int bufoffset = 0; + for (int i = 0; i < sendcount[isend]; i++) { + memcpy(&buf[bufoffset],&sendbuf[nsize*sendindices[m++]],nsize); + bufoffset += nsize; + } + MPI_Send(buf,sendsize[isend],MPI_BYTE,sendproc[isend],0,comm); + } + + // copy self data directly from sendbuf to recvbuf + + if (self) + for (int i = 0; i < sendcount[nsend]; i++) { + memcpy(&recvbuf[recvoffset],&sendbuf[nsize*sendindices[m++]],nsize); + recvoffset += nsize; + } + + // free send buffer + + memory->sfree(buf); + + // wait on all incoming messages + + if (nrecv) MPI_Waitall(nrecv,request,status); +} + +/* ---------------------------------------------------------------------- + sendbuf = data to send + recvbuf = buffer to recv all data into + requires nsendmax,recvsize,sendsize,sendoffset,sendsizedatum + be setup by size(int *, int *) +------------------------------------------------------------------------- */ + +void Irregular::exchange_varying(char *sendbuf, char *recvbuf) +{ + // post all receives + + int recvoffset = 0; + for (int irecv = 0; irecv < nrecv; irecv++) { + MPI_Irecv(&recvbuf[recvoffset],recvsize[irecv],MPI_BYTE, + recvproc[irecv],0,comm,&request[irecv]); + recvoffset += recvsize[irecv]; + } + + // malloc buf for largest send + + char *buf = (char *) memory->smalloc(nsendmax,"buf"); + + // barrier to insure receives are posted + + MPI_Barrier(comm); + + // send each message, packing buf with needed datums + + int index; + int m = 0; + for (int isend = 0; isend < nsend; isend++) { + int bufoffset = 0; + for (int i = 0; i < sendcount[isend]; i++) { + index = sendindices[m++]; + memcpy(&buf[bufoffset],&sendbuf[sendoffset[index]],sendsizedatum[index]); + bufoffset += sendsizedatum[index]; + } + MPI_Send(buf,sendsize[isend],MPI_BYTE,sendproc[isend],0,comm); + } + + // copy self data directly from sendbuf to recvbuf + + if (self) + for (int i = 0; i < sendcount[nsend]; i++) { + index = sendindices[m++]; + memcpy(&recvbuf[recvoffset],&sendbuf[sendoffset[index]], + sendsizedatum[index]); + recvoffset += sendsizedatum[index]; + } + + // free send buffer + + memory->sfree(buf); + + // wait on all incoming messages + + if (nrecv) MPI_Waitall(nrecv,request,status); +} + +/* ---------------------------------------------------------------------- */ + +void Irregular::init() +{ + sendoffsetflag = 0; + sendproc = sendcount = sendsize = sendindices = NULL; + sendoffset = NULL; + recvproc = recvcount = recvsize = NULL; + request = NULL; + status = NULL; +} + +/* ---------------------------------------------------------------------- */ + +void Irregular::deallocate() +{ + delete [] sendproc; + delete [] sendcount; + delete [] sendsize; + memory->sfree(sendindices); + if (sendoffsetflag) memory->sfree(sendoffset); + + delete [] recvproc; + delete [] recvcount; + delete [] recvsize; + + delete [] request; + delete [] status; +} diff --git a/couple/library/irregular.h b/couple/library/irregular.h new file mode 100644 index 0000000000..ea649533ae --- /dev/null +++ b/couple/library/irregular.h @@ -0,0 +1,58 @@ +#ifndef IRREGULAR_H +#define IRREGULAR_H + +#include "mpi.h" + +class Irregular { + public: + Irregular(MPI_Comm); + ~Irregular(); + + void pattern(int, int *); + int size(int); + int size(int *, int *, int *); + void exchange(char *, char *); + + private: + int me,nprocs; + + int patternflag; // UNSET,SET + int sizestyle; // NONE,SAME,VARYING + + int self; // 0 = no data to copy to self, 1 = yes + + int ndatumsend; // # of datums to send w/ self + int ndatumrecv; // # of datums to recv w/ self + int nbytesrecv; // total bytes in received data w/ self + int nsend; // # of messages to send w/out self + int nrecv; // # of messages to recv w/out self + int nsendmax; // # of bytes in largest send message, w/out self + + int *sendproc; // list of procs to send to w/out self + int *sendcount; // # of datums to send to each proc w/ self + int *sendsize; // # of bytes to send to each proc w/ self + int *sendindices; // indices of datums to send to each proc w/ self + + int nsize; // size of every datum in bytes (SAME) + int *sendsizedatum; // bytes in each datum to send w/ self (VARYING) + int *sendoffset; // byte offset to where each datum starts w/ self + int sendoffsetflag; // 1 if allocated sendoffset, 0 if passed in + + int *recvproc; // list of procs to recv from w/out self + int *recvcount; // # of datums to recv from each proc w/out self + int *recvsize; // # of bytes to recv from each proc w/out self + + MPI_Request *request; // MPI requests for posted recvs + MPI_Status *status; // MPI statuses for Waitall + MPI_Comm comm; // MPI communicator for all communication + + class Memory *memory; + class Error *error; + + void exchange_same(char *, char *); + void exchange_varying(char *, char *); + void init(); + void deallocate(); +}; + +#endif diff --git a/couple/library/lammps_data_write.cpp b/couple/library/lammps_data_write.cpp new file mode 100644 index 0000000000..7490537ee5 --- /dev/null +++ b/couple/library/lammps_data_write.cpp @@ -0,0 +1,249 @@ +#include "stdlib.h" +#include "string.h" +#include "lammps_data_write.h" +#include "memory.h" +#include "error.h" + +#define DELTA 4; + +enum{INT,DOUBLE,DOUBLE2}; + +/* ---------------------------------------------------------------------- */ + +LAMMPSDataWrite::LAMMPSDataWrite(MPI_Comm caller_comm) : Send2One(caller_comm) +{ + outfile = NULL; + + nheader = maxheader = 0; + format = NULL; + headtype = NULL; + ihead = NULL; + dhead = NULL; + ddhead = NULL; + + nper = maxper = 0; + atomtype = NULL; + ivec = NULL; + dvec = NULL; + stride = NULL; +} + +/* ---------------------------------------------------------------------- */ + +LAMMPSDataWrite::~LAMMPSDataWrite() +{ + delete [] outfile; + + for (int i = 0; i < nheader; i++) delete [] format[i]; + memory->sfree(format); + memory->sfree(headtype); + memory->sfree(ihead); + memory->sfree(dhead); + memory->destroy_2d_double_array(ddhead); + + memory->sfree(atomtype); + memory->sfree(ivec); + memory->sfree(dvec); + memory->sfree(stride); +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::pre() +{ + if (me == 0) { + fp = fopen(outfile,"w"); + if (fp == NULL) + error->one("Could not open data_write output file"); + } + + if (me == 0) { + fprintf(fp,"%s","LAMMPS data file\n\n"); + for (int i = 0; i < nheader; i++) + if (headtype[i] == INT) fprintf(fp,format[i],ihead[i]); + else if (headtype[i] == DOUBLE) fprintf(fp,format[i],dhead[i]); + else if (headtype[i] == DOUBLE2) fprintf(fp,format[i], + ddhead[i][0],ddhead[i][1]); + fprintf(fp,"\nAtoms\n\n"); + } +} + +/* ---------------------------------------------------------------------- */ + +int LAMMPSDataWrite::size() +{ + return nper*nlocal*sizeof(double); +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::pack(char *cbuf) +{ + int i,j; + + double *dbuf = (double *) cbuf; + + int m = 0; + for (i = 0; i < nlocal; i++) + for (j = 0; j < nper; j++) { + if (i == 0) { + if (atomtype[j] == 0) + printf("AT %d %d %p %d\n", + atomtype[j],stride[j],ivec[j],ivec[j][0]); + else + printf("AT %d %d %p %g\n", + atomtype[j],stride[j],dvec[j],dvec[j][0]); + } + if (atomtype[j] == INT) dbuf[m++] = ivec[j][i*stride[j]]; + else if (atomtype[j] == DOUBLE) dbuf[m++] = dvec[j][i*stride[j]]; + } +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::process(int nbuf, char *cbuf) +{ + int i,j; + + double *dbuf = (double *) cbuf; + int n = nbuf/nper/sizeof(double); + + int m = 0; + for (i = 0; i < n; i++) { + for (j = 0; j < nper; j++) { + double dvalue = dbuf[m++]; + if (atomtype[j] == INT) fprintf(fp,"%d ",static_cast (dvalue)); + else if (atomtype[j] == DOUBLE) fprintf(fp,"%g ",dvalue); + } + fprintf(fp,"\n"); + } +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::post() +{ + if (me == 0) fclose(fp); +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::file(char *fname) +{ + int n = strlen(fname) + 1; + outfile = new char[n]; + strcpy(outfile,fname); +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::header(char *str, int ivalue) +{ + if (nheader == maxheader) grow_header(); + int n = strlen(str) + 2; + format[nheader] = new char[n]; + strcpy(format[nheader],str); + format[nheader][n-2] = '\n'; + format[nheader][n-1] = '\0'; + headtype[nheader] = INT; + ihead[nheader] = ivalue; + nheader++; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::header(char *str, double dvalue) +{ + if (nheader == maxheader) grow_header(); + int n = strlen(str) + 2; + format[nheader] = new char[n]; + strcpy(format[nheader],str); + format[nheader][n-2] = '\n'; + format[nheader][n-1] = '\0'; + headtype[nheader] = DOUBLE; + dhead[nheader] = dvalue; + nheader++; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::header(char *str, double dvalue1, double dvalue2) +{ + if (nheader == maxheader) grow_header(); + int n = strlen(str) + 2; + format[nheader] = new char[n]; + strcpy(format[nheader],str); + format[nheader][n-2] = '\n'; + format[nheader][n-1] = '\0'; + headtype[nheader] = DOUBLE2; + ddhead[nheader][0] = dvalue1; + ddhead[nheader][1] = dvalue2; + nheader++; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::atoms(int n) +{ + nlocal = n; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::atoms(int *vec) +{ + if (nper == maxper) grow_peratom(); + atomtype[nper] = INT; + ivec[nper] = vec; + stride[nper] = 1; + nper++; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::atoms(double *vec) +{ + if (nper == maxper) grow_peratom(); + atomtype[nper] = DOUBLE; + dvec[nper] = vec; + stride[nper] = 1; + nper++; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::atoms(int n, double **vec) +{ + if (nper+n >= maxper) grow_peratom(); + for (int i = 0; i < n; i++) { + atomtype[nper] = DOUBLE; + dvec[nper] = &vec[0][i]; + stride[nper] = n; + nper++; + } +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::grow_header() +{ + int n = maxheader + DELTA; + format = (char **) memory->srealloc(format,n*sizeof(char *),"ldw:format"); + headtype = (int *) memory->srealloc(headtype,n*sizeof(int),"ldw:headtype"); + ihead = (int *) memory->srealloc(ihead,n*sizeof(int),"ldw:ihead"); + dhead = (double *) memory->srealloc(dhead,n*sizeof(double),"ldw:dhead"); + ddhead = memory->grow_2d_double_array(ddhead,n,2,"ldw:ddhead"); + maxheader = n; +} + +/* ---------------------------------------------------------------------- */ + +void LAMMPSDataWrite::grow_peratom() +{ + int n = maxper + DELTA; + atomtype = (int *) memory->srealloc(atomtype,n*sizeof(int *),"ldw:atomtype"); + ivec = (int **) memory->srealloc(ivec,n*sizeof(int *),"ldw:ihead"); + dvec = (double **) memory->srealloc(dvec,n*sizeof(double *),"ldw:dhead"); + stride = (int *) memory->srealloc(stride,n*sizeof(int *),"ldw:stride"); + maxper = n; +} diff --git a/couple/library/lammps_data_write.h b/couple/library/lammps_data_write.h new file mode 100644 index 0000000000..5eda0f4044 --- /dev/null +++ b/couple/library/lammps_data_write.h @@ -0,0 +1,48 @@ +#ifndef LAMMPS_DATA_WRITE_H +#define LAMMPS_DATA_WRITE_H + +#include "send2one.h" +#include "stdio.h" + +class LAMMPSDataWrite : public Send2One { + public: + LAMMPSDataWrite(MPI_Comm); + ~LAMMPSDataWrite(); + + void pre(); + int size(); + void pack(char *); + void process(int, char *); + void post(); + + void file(char *); + void header(char *, int); + void header(char *, double); + void header(char *, double, double); + void atoms(int); + void atoms(int *); + void atoms(double *); + void atoms(int, double **); + + private: + char *outfile; + int nlocal; + FILE *fp; + + int nheader,maxheader; + char **format; + int *headtype,*ihead; + double *dhead; + double **ddhead; + + int nper,maxper; + int *atomtype; + int **ivec; + double **dvec; + int *stride; + + void grow_header(); + void grow_peratom(); +}; + +#endif diff --git a/couple/library/many2many.cpp b/couple/library/many2many.cpp new file mode 100644 index 0000000000..69ed866013 --- /dev/null +++ b/couple/library/many2many.cpp @@ -0,0 +1,316 @@ +#include "mpi.h" +#include "many2many.h" +#include "irregular.h" +#include "memory.h" +#include "error.h" + +#include + +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +Many2Many::Many2Many(MPI_Comm caller) +{ + comm = caller; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + memory = new Memory(comm); + error = new Error(comm); + + src_own = dest_own = NULL; + src_off = dest_off = NULL; + src_iwork = dest_iwork = NULL; + src_dwork = dest_dwork = NULL; + irregular = NULL; +} + +/* ---------------------------------------------------------------------- */ + +Many2Many::~Many2Many() +{ + delete memory; + delete error; + deallocate(); +} + +/* ---------------------------------------------------------------------- + create a many2many pattern, deallocating any previous pattern + each proc will contribute nsrc items with IDs listed in id_src + each proc will receive ndest items with IDs listed in id_dest + only sets up communication via rendezvous algorithm and Irregular class + if id_src does not match id_dest on all procs +------------------------------------------------------------------------- */ + +void Many2Many::setup(int nsrc, int *id_src, int ndest, int *id_dest) +{ + int i,j,isrc,idest,nsend,nrecv; + int *proclist,*work; + std::map hash; + std::map::iterator loc; + + // free any previous many2many info + + deallocate(); + + // allocate on-proc and off-proc index lists + + src_own = + (int *) memory->smalloc(nsrc*sizeof(int),"many2many:src_own"); + dest_own = + (int *) memory->smalloc(ndest*sizeof(int),"many2many:dest_own"); + src_off = + (int *) memory->smalloc(nsrc*sizeof(int),"many2many:src_off"); + dest_off = + (int *) memory->smalloc(ndest*sizeof(int),"many2many:dest_off"); + + // store destination IDs in hash + + for (int i = 0; i < ndest; i++) + hash.insert(std::pair (id_dest[i],i)); + + // src_own, dest_own = list of IDs in both src and dest + // nsrc_off = list of IDs in src owned by other procs + + nown = nsrc_off = 0; + nsrc_off = 0; + for (i = 0; i < nsrc; i++) { + loc = hash.find(id_src[i]); + if (loc != hash.end()) { + src_own[nown] = i; + dest_own[nown] = loc->second; + nown++; + } else src_off[nsrc_off++] = i; + } + + // all done if all procs have one-to-one mapping of src and dest IDs + // else figure out irregular comm needed + + int flag = 0; + if (nown == nsrc && nown == ndest) flag = 1; + int flagall; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MIN,comm); + if (flagall) return; + + // ndest_off = list of IDs in dest owned by other procs + + work = (int *) memory->smalloc(ndest*sizeof(int),"many2many:work"); + + for (i = 0; i < ndest; i++) work[i] = 0; + for (i = 0; i < nown; i++) work[dest_own[i]] = 1; + + ndest_off = 0; + for (i = 0; i < ndest; i++) + if (work[i] == 0) dest_off[ndest_off++] = i; + + memory->sfree(work); + + // realloc off-proc arrays to smaller size + + src_off = (int *) + memory->srealloc(src_off,nsrc_off*sizeof(int),"many2many:src_off"); + dest_off = (int *) + memory->srealloc(dest_off,ndest_off*sizeof(int),"many2many:dest_off"); + + // send off-proc src and dest Datums to 3rd-party proc via irregular comm + // proc = ID % nprocs + + nsend = nsrc_off + ndest_off; + proclist = new int[nsend]; + Datum1 *send1 = new Datum1[nsend]; + + for (i = 0; i < nsrc_off; i++) { + proclist[i] = id_src[src_off[i]] % nprocs; + send1[i].id = id_src[src_off[i]]; + send1[i].proc = me; + send1[i].index = src_off[i]; + } + for (i = 0, j = nsrc_off; i < ndest_off; i++, j++) { + proclist[j] = id_dest[dest_off[i]] % nprocs; + send1[j].id = -id_dest[dest_off[i]]; + send1[j].proc = me; + send1[j].index = dest_off[i]; + } + + irregular = new Irregular(comm); + irregular->pattern(nsend,proclist); + nrecv = irregular->size(sizeof(Datum1)) / sizeof(Datum1); + Datum1 *recv1 = new Datum1[nrecv]; + irregular->exchange((char *) send1, (char *) recv1); + delete irregular; + delete [] proclist; + + // as 3rd-party proc, now have matching pairs of off-proc IDs + // store src IDs (which are positive) in hash + // loop over dest IDs (which are negative) to find matches + // send match info back to src procs via a 2nd irregular comm + + nsend = nrecv/2; + proclist = new int[nsend]; + Datum2 *send2 = new Datum2[nsend]; + nsend = 0; + + hash.clear(); + for (isrc = 0; isrc < nrecv; isrc++) + if (recv1[isrc].id > 0) + hash.insert(std::pair (recv1[isrc].id,isrc)); + + for (idest = 0; idest < nrecv; idest++) + if (recv1[idest].id < 0) { + loc = hash.find(-recv1[idest].id); + if (loc != hash.end()) { + isrc = loc->second; + proclist[nsend] = recv1[isrc].proc; + send2[nsend].slocal = recv1[isrc].index; + send2[nsend].dlocal = recv1[idest].index; + send2[nsend].dproc = recv1[idest].proc; + nsend++; + } else error->one("Did not receive matching src/dest ID"); + } + + irregular = new Irregular(comm); + irregular->pattern(nsend,proclist); + nrecv = irregular->size(sizeof(Datum2)) / sizeof(Datum2); + Datum2 *recv2 = new Datum2[nrecv]; + irregular->exchange((char *) send2, (char *) recv2); + delete irregular; + delete [] proclist; + + // use list of received src->dest Datums to build final irregular commm + // irregular comm will communicate off-proc info from src to dest directly + // work = local indices of dest IDs to send initially + + nsend = nrecv; + proclist = new int[nsend]; + work = new int[nsend]; + + for (i = 0; i < nrecv; i++) { + src_off[i] = recv2[i].slocal; + work[i] = recv2[i].dlocal; + proclist[i] = recv2[i].dproc; + } + + irregular = new Irregular(comm); + irregular->pattern(nsend,proclist); + + // send receiver's local indices + // receiver stores them as indirection list in dest_off + + nrecv = irregular->size(sizeof(int)) / sizeof(int); + irregular->exchange((char *) work, (char *) dest_off); + + delete [] proclist; + delete [] work; + + // create work arrays for data exchange of int/double data + + src_iwork = + (int *) memory->smalloc(nsrc_off*sizeof(int),"many2many:src_iwork"); + dest_iwork = + (int *) memory->smalloc(ndest_off*sizeof(int),"many2many:dest_iwork"); + src_dwork = + (double *) memory->smalloc(nsrc_off*sizeof(double),"many2many:src_dwork"); + dest_dwork = + (double *) memory->smalloc(ndest_off*sizeof(double), + "many2many:dest_dwork"); + + // clean up + + delete [] send1; + delete [] recv1; + delete [] send2; + delete [] recv2; + + // debug checks for full coverage of srd/dest - can delete eventually + + work = new int[MAX(nsrc,ndest)]; + + for (i = 0; i < nsrc; i++) work[i] = 0; + for (i = 0; i < nown; i++) work[src_own[i]]++; + for (i = 0; i < nsrc_off; i++) work[src_off[i]]++; + for (i = 0; i < nsrc; i++) + if (work[i] != 1) { + char str[128]; + sprintf(str,"BAD SRC %d: %d %d\n",me,i,work[i]); + error->one(str); + } + + for (i = 0; i < ndest; i++) work[i] = 0; + for (i = 0; i < nown; i++) work[dest_own[i]]++; + for (i = 0; i < ndest_off; i++) work[dest_off[i]]++; + for (i = 0; i < ndest; i++) + if (work[i] != 1) { + char str[128]; + sprintf(str,"BAD DEST %d: %d %d\n",me,i,work[i]); + error->one(str); + } + + delete [] work; +} + +/* ---------------------------------------------------------------------- + transfer one src entity to dest entity, matched by IDs in create() + operates on an int vector +------------------------------------------------------------------------- */ + +void Many2Many::exchange(int *src, int *dest) +{ + int i; + + // copy on-proc info + + for (i = 0; i < nown; i++) + dest[dest_own[i]] = src[src_own[i]]; + + // communicate off-proc info + // user src_off and dest_off to pack/unpack data + + if (irregular) { + int nrecv = irregular->size(sizeof(int)) / sizeof(int); + for (i = 0; i < nsrc_off; i++) src_iwork[i] = src[src_off[i]]; + irregular->exchange((char *) src_iwork, (char *) dest_iwork); + for (i = 0; i < ndest_off; i++) dest[dest_off[i]] = dest_iwork[i]; + } +} + +/* ---------------------------------------------------------------------- + transfer one src entity to dest entity, matched by IDs in create() + operates on a double vector +------------------------------------------------------------------------- */ + +void Many2Many::exchange(double *src, double *dest) +{ + int i; + + // copy on-proc info + + for (int i = 0; i < nown; i++) + dest[dest_own[i]] = src[src_own[i]]; + + // communicate off-proc info + // user src_off and dest_off to pack/unpack data + + if (irregular) { + int nrecv = irregular->size(sizeof(double)) / sizeof(double); + for (i = 0; i < nsrc_off; i++) src_dwork[i] = src[src_off[i]]; + irregular->exchange((char *) src_dwork, (char *) dest_dwork); + for (i = 0; i < ndest_off; i++) dest[dest_off[i]] = dest_dwork[i]; + } +} + +/* ---------------------------------------------------------------------- */ + +void Many2Many::deallocate() +{ + memory->sfree(src_own); + memory->sfree(dest_own); + memory->sfree(src_off); + memory->sfree(dest_off); + memory->sfree(src_iwork); + memory->sfree(dest_iwork); + memory->sfree(src_dwork); + memory->sfree(dest_dwork); + + delete irregular; +} diff --git a/couple/library/many2many.h b/couple/library/many2many.h new file mode 100644 index 0000000000..a0c4a2b4db --- /dev/null +++ b/couple/library/many2many.h @@ -0,0 +1,47 @@ +#ifndef MANY2MANY_H +#define MANY2MANY_H + +#include "mpi.h" + +class Many2Many { + public: + Many2Many(MPI_Comm); + ~Many2Many(); + + void setup(int, int *, int, int *); + void exchange(int *, int *); + void exchange(double *, double *); + + protected: + int me,nprocs; + MPI_Comm comm; + class Memory *memory; + class Error *error; + + int nown; // # of IDs common to src and dest + int nsrc_off,ndest_off; // # of off-processor IDs + + int *src_own,*dest_own; // indices of the owned IDs + int *src_off,*dest_off; // indices of the off-proc IDs + + int *src_iwork,*dest_iwork; // work arrays for comm of ints + double *src_dwork,*dest_dwork; // work arrays for comm of doubles + + class Irregular *irregular; // irregular comm from src->dest + + struct Datum1 { + int id; // src or dest global ID + int proc; // owning proc + int index; // local index on owning proc + }; + + struct Datum2 { + int slocal; // local index of src ID on sending proc + int dlocal; // local index of dest ID on receiving proc + int dproc; // receiving proc + }; + + void deallocate(); +}; + +#endif diff --git a/couple/library/many2one.cpp b/couple/library/many2one.cpp new file mode 100644 index 0000000000..5b5467aa61 --- /dev/null +++ b/couple/library/many2one.cpp @@ -0,0 +1,103 @@ +#include "mpi.h" +#include "stdio.h" +#include "stdlib.h" +#include "many2one.h" +#include "memory.h" + +/* ---------------------------------------------------------------------- */ + +Many2One::Many2One(MPI_Comm caller_comm) +{ + comm = caller_comm; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + memory = new Memory(comm); + + if (me == 0) { + counts = new int[nprocs]; + multicounts = new int[nprocs]; + displs = new int[nprocs]; + multidispls = new int[nprocs]; + } else counts = multicounts = displs = multidispls = NULL; + + idall = NULL; +} + +/* ---------------------------------------------------------------------- */ + +Many2One::~Many2One() +{ + delete memory; + + delete [] counts; + delete [] multicounts; + delete [] displs; + delete [] multidispls; + + memory->sfree(idall); +} + +/* ---------------------------------------------------------------------- */ + +void Many2One::setup(int nsrc_in, int *id, int ndest) +{ + nsrc = nsrc_in; + MPI_Allreduce(&nsrc,&nall,1,MPI_INT,MPI_SUM,comm); + MPI_Gather(&nsrc,1,MPI_INT,counts,1,MPI_INT,0,comm); + + if (me == 0) { + displs[0] = 0; + for (int i = 1; i < nprocs; i++) + displs[i] = displs[i-1] + counts[i-1]; + } + + // gather IDs into idall + + idall = NULL; + if (me == 0) + idall = (int *) memory->smalloc(nall*sizeof(int),"many2one:idall"); + MPI_Gatherv(id,nsrc,MPI_INT,idall,counts,displs,MPI_INT,0,comm); +} + +/* ---------------------------------------------------------------------- */ + +void Many2One::gather(double *src, int n, double *dest) +{ + int i,j,ii,jj,m; + + if (me == 0) + for (int i = 0; i < nprocs; i++) { + multicounts[i] = n*counts[i]; + multidispls[i] = n*displs[i]; + } + + // allgather src into desttmp + + double *desttmp = NULL; + if (me == 0) + desttmp = (double *) memory->smalloc(n*nall*sizeof(double), + "many2one:idsttmp"); + MPI_Gatherv(src,n*nsrc,MPI_DOUBLE,desttmp,multicounts,multidispls, + MPI_DOUBLE,0,comm); + + // use idall to move datums from desttmp to dest + + if (me == 0) { + if (n == 1) + for (i = 0; i < nall; i++) { + j = idall[i] - 1; + dest[j] = desttmp[i]; + } + else + for (i = 0; i < nall; i++) { + j = idall[i] - 1; + ii = n*i; + jj = n*j; + for (m = 0; m < n; m++) + dest[jj++] = desttmp[ii++]; + } + } + + memory->sfree(desttmp); +} diff --git a/couple/library/many2one.h b/couple/library/many2one.h new file mode 100644 index 0000000000..7e487ec473 --- /dev/null +++ b/couple/library/many2one.h @@ -0,0 +1,25 @@ +#ifndef MANY2ONE_H +#define MANY2ONE_H + +#include "mpi.h" + +class Many2One { + public: + Many2One(MPI_Comm); + ~Many2One(); + + void setup(int, int *, int); + void gather(double *, int, double *); + + protected: + int me,nprocs; + MPI_Comm comm; + class Memory *memory; + + int nsrc,nall; + int *counts,*multicounts; + int *displs,*multidispls; + int *idall; +}; + +#endif diff --git a/couple/library/memory.cpp b/couple/library/memory.cpp new file mode 100644 index 0000000000..7377fac744 --- /dev/null +++ b/couple/library/memory.cpp @@ -0,0 +1,120 @@ +#include "mpi.h" +#include "stdlib.h" +#include "stdio.h" +#include "memory.h" +#include "error.h" + +/* ---------------------------------------------------------------------- */ + +Memory::Memory(MPI_Comm comm) +{ + error = new Error(comm); +} + +/* ---------------------------------------------------------------------- */ + +Memory::~Memory() +{ + delete error; +} + +/* ---------------------------------------------------------------------- + safe malloc +------------------------------------------------------------------------- */ + +void *Memory::smalloc(int n, const char *name) +{ + if (n == 0) return NULL; + void *ptr = malloc(n); + if (ptr == NULL) { + char str[128]; + sprintf(str,"Failed to allocate %d bytes for array %s",n,name); + error->one(str); + } + return ptr; +} + +/* ---------------------------------------------------------------------- + safe free +------------------------------------------------------------------------- */ + +void Memory::sfree(void *ptr) +{ + if (ptr == NULL) return; + free(ptr); +} + +/* ---------------------------------------------------------------------- + safe realloc +------------------------------------------------------------------------- */ + +void *Memory::srealloc(void *ptr, int n, const char *name) +{ + if (n == 0) { + sfree(ptr); + return NULL; + } + + ptr = realloc(ptr,n); + if (ptr == NULL) { + char str[128]; + sprintf(str,"Failed to reallocate %d bytes for array %s",n,name); + error->one(str); + } + return ptr; +} + +/* ---------------------------------------------------------------------- + create a 2d double array +------------------------------------------------------------------------- */ + +double **Memory::create_2d_double_array(int n1, int n2, const char *name) + +{ + double *data = (double *) smalloc(n1*n2*sizeof(double),name); + double **array = (double **) smalloc(n1*sizeof(double *),name); + + int n = 0; + for (int i = 0; i < n1; i++) { + array[i] = &data[n]; + n += n2; + } + + return array; +} + +/* ---------------------------------------------------------------------- + grow or shrink 1st dim of a 2d double array + last dim must stay the same + if either dim is 0, return NULL +------------------------------------------------------------------------- */ + +double **Memory::grow_2d_double_array(double **array, + int n1, int n2, const char *name) + +{ + if (array == NULL) return create_2d_double_array(n1,n2,name); + + double *data = (double *) srealloc(array[0],n1*n2*sizeof(double),name); + array = (double **) srealloc(array,n1*sizeof(double *),name); + + int n = 0; + for (int i = 0; i < n1; i++) { + array[i] = &data[n]; + n += n2; + } + + return array; +} + +/* ---------------------------------------------------------------------- + free a 2d double array +------------------------------------------------------------------------- */ + +void Memory::destroy_2d_double_array(double **array) + +{ + if (array == NULL) return; + sfree(array[0]); + sfree(array); +} diff --git a/couple/library/memory.h b/couple/library/memory.h new file mode 100644 index 0000000000..d6efa7789e --- /dev/null +++ b/couple/library/memory.h @@ -0,0 +1,23 @@ +#ifndef MEMORY_H +#define MEMORY_H + +#include "mpi.h" + +class Memory { + public: + Memory(MPI_Comm); + ~Memory(); + + void *smalloc(int n, const char *); + void sfree(void *); + void *srealloc(void *, int n, const char *name); + + double **create_2d_double_array(int, int, const char *); + double **grow_2d_double_array(double **, int, int, const char *); + void destroy_2d_double_array(double **); + + private: + class Error *error; +}; + +#endif diff --git a/couple/library/one2many.cpp b/couple/library/one2many.cpp new file mode 100644 index 0000000000..e26bf7086f --- /dev/null +++ b/couple/library/one2many.cpp @@ -0,0 +1,76 @@ +#include "mpi.h" +#include "one2many.h" +#include "memory.h" + +#include + +/* ---------------------------------------------------------------------- */ + +One2Many::One2Many(MPI_Comm caller_comm) +{ + comm = caller_comm; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + memory = new Memory(comm); + hash = new std::map(); +} + +/* ---------------------------------------------------------------------- */ + +One2Many::~One2Many() +{ + delete memory; + delete hash; +} + +/* ---------------------------------------------------------------------- */ + +void One2Many::setup(int nsrc_in, int ndest, int *id) +{ + nsrc = nsrc_in; + + // store my local IDs in hash + + hash->clear(); + for (int i = 0; i < ndest; i++) + hash->insert(std::pair (id[i],i)); +} + +/* ---------------------------------------------------------------------- */ + +void One2Many::scatter(double *src, int n, double *dest) +{ + int i,j,k,m; + + // allocate src on procs that don't have it + + int flag = 0; + if (src == NULL) { + src = (double *) memory->smalloc(n*nsrc*sizeof(double),"one2many:src"); + flag = 1; + } + + // broadcast src from 0 to other procs + + MPI_Bcast(src,n*nsrc,MPI_DOUBLE,0,comm); + + // each proc loops over entire src + // if I own the global ID, copy src values into dest + + std::map::iterator loc; + for (m = 1; m <= nsrc; m++) { + loc = hash->find(m); + if (loc == hash->end()) continue; + i = n*loc->second; + j = 3*(m-1); + if (n == 1) dest[i] = src[j]; + else + for (k = 0; k < n; k++) + dest[i++] = src[j++]; + } + + // free locally allocated src + + if (flag) memory->sfree(src); +} diff --git a/couple/library/one2many.h b/couple/library/one2many.h new file mode 100644 index 0000000000..9797866112 --- /dev/null +++ b/couple/library/one2many.h @@ -0,0 +1,24 @@ +#ifndef ONE2MANY_H +#define ONE2MANY_H + +#include "mpi.h" + +#include + +class One2Many { + public: + One2Many(MPI_Comm); + ~One2Many(); + + void setup(int, int, int *); + void scatter(double *, int, double *); + + protected: + int me,nprocs; + MPI_Comm comm; + class Memory *memory; + std::map *hash; + int nsrc; +}; + +#endif diff --git a/couple/library/send2one.cpp b/couple/library/send2one.cpp new file mode 100644 index 0000000000..d969b4678d --- /dev/null +++ b/couple/library/send2one.cpp @@ -0,0 +1,83 @@ +#include "mpi.h" +#include "stdlib.h" +#include "stdio.h" +#include "send2one.h" +#include "memory.h" +#include "error.h" + +/* ---------------------------------------------------------------------- */ + +Send2One::Send2One(MPI_Comm caller_comm) +{ + comm = caller_comm; + MPI_Comm_rank(comm,&me); + MPI_Comm_size(comm,&nprocs); + + memory = new Memory(comm); + error = new Error(comm); + + buf = NULL; + maxbuf = 0; +} + +/* ---------------------------------------------------------------------- */ + +Send2One::~Send2One() +{ + delete memory; + delete error; + memory->sfree(buf); +} + +/* ---------------------------------------------------------------------- */ + +void Send2One::execute() +{ + int nme,nmax,nsize,ping; + MPI_Status status; + MPI_Request request; + + // pre-processing before ping loop + + pre(); + + // nme = size of data I contribute, in bytes + // nmax = max size of data on any proc, in bytes + // reallocate buf if necessary + + nme = size(); + MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,comm); + + if (nmax > maxbuf) { + maxbuf = nmax; + memory->sfree(buf); + buf = (char *) memory->smalloc(maxbuf,"foo:buf"); + } + + // pack my data into buf + + pack(buf); + + // proc 0 pings each proc, receives its data + // all other procs wait for ping, send their data to proc 0 + // invoke process() to work with data + + if (me == 0) { + for (int iproc = 0; iproc < nprocs; iproc++) { + if (iproc) { + MPI_Irecv(buf,maxbuf,MPI_CHAR,iproc,0,comm,&request); + MPI_Send(&ping,0,MPI_INT,iproc,0,comm); + MPI_Wait(&request,&status); + MPI_Get_count(&status,MPI_CHAR,&nsize); + } else nsize = nme; + + process(nsize,buf); + } + + } else { + MPI_Recv(&ping,0,MPI_INT,0,0,comm,&status); + MPI_Rsend(buf,nme,MPI_CHAR,0,0,comm); + } + + post(); +} diff --git a/couple/library/send2one.h b/couple/library/send2one.h new file mode 100644 index 0000000000..df25b5fb70 --- /dev/null +++ b/couple/library/send2one.h @@ -0,0 +1,29 @@ +#ifndef SEND2ONE_H +#define SEND2ONE_H + +#include "mpi.h" + +class Send2One { + public: + Send2One(MPI_Comm); + virtual ~Send2One(); + + void execute(); + + protected: + int me,nprocs; + MPI_Comm comm; + class Memory *memory; + class Error *error; + + int maxbuf; + char *buf; + + virtual void pre() = 0; + virtual int size() = 0; + virtual void pack(char *) = 0; + virtual void process(int, char *) = 0; + virtual void post() = 0; +}; + +#endif diff --git a/couple/simple/README b/couple/simple/README new file mode 100644 index 0000000000..507596d77f --- /dev/null +++ b/couple/simple/README @@ -0,0 +1,58 @@ +This directory has a simple C and C++ code that shows how LAMMPS can +be linked to a driver application as a library. The purpose is to +illustrate how another code could perform computations while using +LAMMPS to perform MD, or how an umbrella code or script could call +both LAMMPS and some other code to perform a coupled calculation. + +simple.cpp is the C++ driver +simple.c is the C driver + +The 2 codes do the same thing, so you can compare them to see how to +drive LAMMPS in this manner. The C driver is similar in spirit to +what one could use from a Fortran program or scripting language. + +You can then build either driver code with a compile line something +like this, which includes paths to the LAMMPS library interface, MPI, +and FFTW (assuming you built LAMMPS as a library with its PPPM +solver). + +This builds the C++ driver with the LAMMPS library using a C++ compiler: + +g++ -I/home/sjplimp/lammps/src -c simple.cpp +g++ -L/home/sjplimp/lammps/src simple.o \ + -llmp_g++ -lfftw -lmpich -lpthread -o simpleCC + +This builds the C driver with the LAMMPS library using a C compiler: + +gcc -I/home/sjplimp/lammps/src -c simple.c +gcc -L/home/sjplimp/lammps/src simple.o \ + -llmp_g++ -lfftw -lmpich -lpthread -lstdc++ -o simpleC + +You then run simpleCC or simpleC on a parallel machine on some number +of processors Q with 2 arguments: + +mpirun -np Q simpleCC P in.lj + +P is the number of procs you want LAMMPS to run on (must be <= Q) and +in.lj is a LAMMPS input script. + +The driver will launch LAMMPS on P procs, read the input script a line +at a time, and pass each command line to LAMMPS. The final line of +the script is a "run" command, so LAMMPS will run the problem. + +The driver then requests all the atom coordinates from LAMMPS, moves +one of the atoms a small amount "epsilon", passes the coordinates back +to LAMMPS, and runs LAMMPS again. If you look at the output, you +should see a small energy change between runs, due to the moved atom. + +The C driver is calling C-style routines in the src/library.cpp file +of LAMMPS. You could add any functions you wish to this file to +manipulate LAMMPS data however you wish. + +The C++ driver does the same thing, except that it instantiates LAMMPS +as an object first. Some of the functions in src/library.cpp can be +invoked directly as methods within appropriate LAMMPS classes, which +is what the driver does. Any public LAMMPS class method could be +called from the driver this way. However the get/put functions are +only implemented in src/library.cpp, so the C++ driver calls them as +C-style functions. diff --git a/couple/simple/in.lj b/couple/simple/in.lj new file mode 100644 index 0000000000..98de6750b4 --- /dev/null +++ b/couple/simple/in.lj @@ -0,0 +1,24 @@ +# 3d Lennard-Jones melt + +units lj +atom_style atomic +atom_modify map array + +lattice fcc 0.8442 +region box block 0 4 0 4 0 4 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 20 check no + +fix 1 all nve + +run 10 + diff --git a/couple/simple/simple.c b/couple/simple/simple.c new file mode 100644 index 0000000000..5b80815160 --- /dev/null +++ b/couple/simple/simple.c @@ -0,0 +1,116 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +/* c_driver = simple example of how an umbrella program + can invoke LAMMPS as a library on some subset of procs + Syntax: c_driver P in.lammps + P = # of procs to run LAMMPS on + must be <= # of procs the driver code itself runs on + in.lammps = LAMMPS input script + See README for compilation instructions */ + +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "mpi.h" +#include "library.h" /* this is a LAMMPS include file */ + +int main(int narg, char **arg) +{ + /* setup MPI and various communicators + driver runs on all procs in MPI_COMM_WORLD + comm_lammps only has 1st P procs (could be all or any subset) */ + + MPI_Init(&narg,&arg); + + if (narg != 3) { + printf("Syntax: c_driver P in.lammps\n"); + exit(1); + } + + int me,nprocs; + MPI_Comm_rank(MPI_COMM_WORLD,&me); + MPI_Comm_size(MPI_COMM_WORLD,&nprocs); + + int nprocs_lammps = atoi(arg[1]); + if (nprocs_lammps > nprocs) { + if (me == 0) + printf("ERROR: LAMMPS cannot use more procs than available\n"); + MPI_Abort(MPI_COMM_WORLD,1); + } + + int lammps; + if (me < nprocs_lammps) lammps = 1; + else lammps = MPI_UNDEFINED; + MPI_Comm comm_lammps; + MPI_Comm_split(MPI_COMM_WORLD,lammps,0,&comm_lammps); + + /* open LAMMPS input script */ + + FILE *fp; + if (me == 0) { + fp = fopen(arg[2],"r"); + if (fp == NULL) { + printf("ERROR: Could not open LAMMPS input script\n"); + MPI_Abort(MPI_COMM_WORLD,1); + } + } + + /* run the input script thru LAMMPS one line at a time until end-of-file + driver proc 0 reads a line, Bcasts it to all procs + (could just send it to proc 0 of comm_lammps and let it Bcast) + all LAMMPS procs call lammps_command() on the line */ + + void *ptr; + if (lammps == 1) lammps_open(0,NULL,comm_lammps,&ptr); + + int n; + char line[1024]; + while (1) { + if (me == 0) { + if (fgets(line,1024,fp) == NULL) n = 0; + else n = strlen(line) + 1; + if (n == 0) fclose(fp); + } + MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD); + if (n == 0) break; + MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD); + if (lammps == 1) lammps_command(ptr,line); + } + + /* run 10 more steps + get coords from LAMMPS + change coords of 1st atom + put coords back into LAMMPS + run a single step with changed coords */ + + if (lammps == 1) { + lammps_command(ptr,"run 10"); + + int natoms = lammps_get_natoms(ptr); + double *x = (double *) malloc(3*natoms*sizeof(double)); + lammps_get_coords(ptr,x); + double epsilon = 0.1; + x[0] += epsilon; + lammps_put_coords(ptr,x); + free(x); + + lammps_command(ptr,"run 1"); + } + + if (lammps == 1) lammps_close(ptr); + + /* close down MPI */ + + MPI_Finalize(); +} diff --git a/couple/simple/simple.cpp b/couple/simple/simple.cpp new file mode 100644 index 0000000000..ff78e49657 --- /dev/null +++ b/couple/simple/simple.cpp @@ -0,0 +1,121 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +// c++_driver = simple example of how an umbrella program +// can invoke LAMMPS as a library on some subset of procs +// Syntax: c++_driver P in.lammps +// P = # of procs to run LAMMPS on +// must be <= # of procs the driver code itself runs on +// in.lammps = LAMMPS input script +// See README for compilation instructions + +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "mpi.h" +#include "lammps.h" // these are LAMMPS include files +#include "input.h" +#include "atom.h" +#include "library.h" + +using namespace LAMMPS_NS; + +int main(int narg, char **arg) +{ + // setup MPI and various communicators + // driver runs on all procs in MPI_COMM_WORLD + // comm_lammps only has 1st P procs (could be all or any subset) + + MPI_Init(&narg,&arg); + + if (narg != 3) { + printf("Syntax: c++_driver P in.lammps\n"); + exit(1); + } + + int me,nprocs; + MPI_Comm_rank(MPI_COMM_WORLD,&me); + MPI_Comm_size(MPI_COMM_WORLD,&nprocs); + + int nprocs_lammps = atoi(arg[1]); + if (nprocs_lammps > nprocs) { + if (me == 0) + printf("ERROR: LAMMPS cannot use more procs than available\n"); + MPI_Abort(MPI_COMM_WORLD,1); + } + + int lammps; + if (me < nprocs_lammps) lammps = 1; + else lammps = MPI_UNDEFINED; + MPI_Comm comm_lammps; + MPI_Comm_split(MPI_COMM_WORLD,lammps,0,&comm_lammps); + + // open LAMMPS input script + + FILE *fp; + if (me == 0) { + fp = fopen(arg[2],"r"); + if (fp == NULL) { + printf("ERROR: Could not open LAMMPS input script\n"); + MPI_Abort(MPI_COMM_WORLD,1); + } + } + + // run the input script thru LAMMPS one line at a time until end-of-file + // driver proc 0 reads a line, Bcasts it to all procs + // (could just send it to proc 0 of comm_lammps and let it Bcast) + // all LAMMPS procs call input->one() on the line + + LAMMPS *lmp; + if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps); + + int n; + char line[1024]; + while (1) { + if (me == 0) { + if (fgets(line,1024,fp) == NULL) n = 0; + else n = strlen(line) + 1; + if (n == 0) fclose(fp); + } + MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD); + if (n == 0) break; + MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD); + if (lammps == 1) lmp->input->one(line); + } + + // run 10 more steps + // get coords from LAMMPS + // change coords of 1st atom + // put coords back into LAMMPS + // run a single step with changed coords + + if (lammps == 1) { + lmp->input->one("run 10"); + + int natoms = static_cast (lmp->atom->natoms); + double *x = new double[3*natoms]; + lammps_get_coords(lmp,x); // no LAMMPS class function for this + double epsilon = 0.1; + x[0] += epsilon; + lammps_put_coords(lmp,x); // no LAMMPS class function for this + delete [] x; + + lmp->input->one("run 1"); + } + + if (lammps == 1) delete lmp; + + // close down MPI + + MPI_Finalize(); +} diff --git a/examples/README b/examples/README index 6212802230..0af97f5d38 100644 --- a/examples/README +++ b/examples/README @@ -62,10 +62,6 @@ You can visualize the dump file as follows: ------------------------------------------ -There is also a COUPLE directory which has examples of how to link to -LAMMPS as a library and drive it from a wrapper program. See the -README file for more info. - There is also an ELASTIC directory with an example script for computing elastic constants, using a zero temperature Si example. See the in.elastic file for more info. diff --git a/src/Makefile.lib b/src/Makefile.lib index 8556846824..6d6cfec7ee 100644 --- a/src/Makefile.lib +++ b/src/Makefile.lib @@ -7,9 +7,9 @@ SHELL = /bin/sh ROOT = lmp EXE = lib$(ROOT)_$@.a -SRC = angle.cpp angle_charmm.cpp angle_cosine.cpp angle_cosine_delta.cpp angle_cosine_squared.cpp angle_harmonic.cpp angle_hybrid.cpp angle_table.cpp atom.cpp atom_vec.cpp atom_vec_angle.cpp atom_vec_atomic.cpp atom_vec_bond.cpp atom_vec_charge.cpp atom_vec_full.cpp atom_vec_hybrid.cpp atom_vec_molecular.cpp bond.cpp bond_fene.cpp bond_fene_expand.cpp bond_harmonic.cpp bond_hybrid.cpp bond_morse.cpp bond_nonlinear.cpp bond_quartic.cpp bond_table.cpp change_box.cpp comm.cpp compute.cpp compute_angle_local.cpp compute_bond_local.cpp compute_centro_atom.cpp compute_cna_atom.cpp compute_com.cpp compute_com_molecule.cpp compute_coord_atom.cpp compute_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_sphere.cpp compute_group_group.cpp compute_gyration.cpp compute_gyration_molecule.cpp compute_heat_flux.cpp compute_improper_local.cpp compute_ke.cpp compute_ke_atom.cpp compute_msd.cpp compute_msd_molecule.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_pressure.cpp compute_property_atom.cpp compute_property_local.cpp compute_property_molecule.cpp compute_rdf.cpp compute_reduce.cpp compute_reduce_region.cpp compute_stress_atom.cpp compute_temp.cpp compute_temp_com.cpp compute_temp_deform.cpp compute_temp_partial.cpp compute_temp_profile.cpp compute_temp_ramp.cpp compute_temp_region.cpp compute_temp_sphere.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.cpp dihedral_harmonic.cpp dihedral_helix.cpp dihedral_hybrid.cpp dihedral_multi_harmonic.cpp dihedral_opls.cpp displace_atoms.cpp displace_box.cpp domain.cpp dump.cpp dump_atom.cpp dump_cfg.cpp dump_custom.cpp dump_dcd.cpp dump_local.cpp dump_xyz.cpp error.cpp ewald.cpp fft3d.cpp fft3d_wrap.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.cpp fix_ave_atom.cpp fix_ave_correlate.cpp fix_ave_histo.cpp fix_ave_spatial.cpp fix_ave_time.cpp fix_aveforce.cpp fix_bond_break.cpp fix_bond_create.cpp fix_bond_swap.cpp fix_box_relax.cpp fix_deform.cpp fix_deposit.cpp fix_drag.cpp fix_dt_reset.cpp fix_efield.cpp fix_enforce2d.cpp fix_evaporate.cpp fix_gravity.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_nh.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_sphere.cpp fix_npt.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_limit.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nvt.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_orient_fcc.cpp fix_planeforce.cpp fix_press_berendsen.cpp fix_print.cpp fix_qeq_comb.cpp fix_recenter.cpp fix_respa.cpp fix_rigid.cpp fix_rigid_nve.cpp fix_rigid_nvt.cpp fix_setforce.cpp fix_shake.cpp fix_shear_history.cpp fix_spring.cpp fix_spring_rg.cpp fix_spring_self.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_harmonic.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_reflect.cpp fix_wall_region.cpp force.cpp group.cpp improper.cpp improper_cvff.cpp improper_harmonic.cpp improper_hybrid.cpp input.cpp integrate.cpp kspace.cpp lammps.cpp lattice.cpp library.cpp memory.cpp min.cpp min_cg.cpp min_hftn.cpp min_linesearch.cpp min_sd.cpp minimize.cpp modify.cpp neigh_bond.cpp neigh_derive.cpp neigh_full.cpp neigh_gran.cpp neigh_half_bin.cpp neigh_half_multi.cpp neigh_half_nsq.cpp neigh_list.cpp neigh_request.cpp neigh_respa.cpp neigh_stencil.cpp neighbor.cpp output.cpp pack.cpp pair.cpp pair_airebo.cpp pair_born_coul_long.cpp pair_buck.cpp pair_buck_coul_cut.cpp pair_buck_coul_long.cpp pair_comb.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_long.cpp pair_dpd.cpp pair_dpd_tstat.cpp pair_eam.cpp pair_eam_alloy.cpp pair_eam_fs.cpp pair_eim.cpp pair_hybrid.cpp pair_hybrid_overlay.cpp pair_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.cpp pair_lj_charmm_coul_long.cpp pair_lj_cut.cpp pair_lj_cut_coul_cut.cpp pair_lj_cut_coul_debye.cpp pair_lj_cut_coul_long.cpp pair_lj_cut_coul_long_tip4p.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_smooth.cpp pair_morse.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_zbl.cpp pair_yukawa.cpp pppm.cpp pppm_tip4p.cpp random_mars.cpp random_park.cpp read_data.cpp read_restart.cpp region.cpp region_block.cpp region_cone.cpp region_cylinder.cpp region_intersect.cpp region_plane.cpp region_prism.cpp region_sphere.cpp region_union.cpp remap.cpp remap_wrap.cpp replicate.cpp respa.cpp run.cpp set.cpp shell.cpp special.cpp temper.cpp thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp write_restart.cpp +SRC = angle.cpp angle_charmm.cpp angle_class2.cpp angle_cosine.cpp angle_cosine_delta.cpp angle_cosine_squared.cpp angle_harmonic.cpp angle_hybrid.cpp angle_table.cpp atom.cpp atom_vec.cpp atom_vec_angle.cpp atom_vec_atomic.cpp atom_vec_bond.cpp atom_vec_charge.cpp atom_vec_colloid.cpp atom_vec_dipole.cpp atom_vec_ellipsoid.cpp atom_vec_full.cpp atom_vec_granular.cpp atom_vec_hybrid.cpp atom_vec_molecular.cpp atom_vec_peri.cpp bond.cpp bond_class2.cpp bond_fene.cpp bond_fene_expand.cpp bond_harmonic.cpp bond_hybrid.cpp bond_morse.cpp bond_nonlinear.cpp bond_quartic.cpp bond_table.cpp change_box.cpp comm.cpp compute.cpp compute_angle_local.cpp compute_bond_local.cpp compute_centro_atom.cpp compute_cna_atom.cpp compute_com.cpp compute_com_molecule.cpp compute_coord_atom.cpp compute_damage_atom.cpp compute_dihedral_local.cpp compute_displace_atom.cpp compute_erotate_asphere.cpp compute_erotate_sphere.cpp compute_event_displace.cpp compute_group_group.cpp compute_gyration.cpp compute_gyration_molecule.cpp compute_heat_flux.cpp compute_improper_local.cpp compute_ke.cpp compute_ke_atom.cpp compute_msd.cpp compute_msd_molecule.cpp compute_pair_local.cpp compute_pe.cpp compute_pe_atom.cpp compute_pressure.cpp compute_property_atom.cpp compute_property_local.cpp compute_property_molecule.cpp compute_rdf.cpp compute_reduce.cpp compute_reduce_region.cpp compute_stress_atom.cpp compute_temp.cpp compute_temp_asphere.cpp compute_temp_com.cpp compute_temp_deform.cpp compute_temp_partial.cpp compute_temp_profile.cpp compute_temp_ramp.cpp compute_temp_region.cpp compute_temp_sphere.cpp create_atoms.cpp create_box.cpp delete_atoms.cpp delete_bonds.cpp dihedral.cpp dihedral_charmm.cpp dihedral_class2.cpp dihedral_harmonic.cpp dihedral_helix.cpp dihedral_hybrid.cpp dihedral_multi_harmonic.cpp dihedral_opls.cpp displace_atoms.cpp displace_box.cpp domain.cpp dump.cpp dump_atom.cpp dump_cfg.cpp dump_custom.cpp dump_dcd.cpp dump_local.cpp dump_xtc.cpp dump_xyz.cpp error.cpp ewald.cpp fft3d.cpp fft3d_wrap.cpp finish.cpp fix.cpp fix_adapt.cpp fix_addforce.cpp fix_ave_atom.cpp fix_ave_correlate.cpp fix_ave_histo.cpp fix_ave_spatial.cpp fix_ave_time.cpp fix_aveforce.cpp fix_bond_break.cpp fix_bond_create.cpp fix_bond_swap.cpp fix_box_relax.cpp fix_deform.cpp fix_deposit.cpp fix_drag.cpp fix_dt_reset.cpp fix_efield.cpp fix_enforce2d.cpp fix_evaporate.cpp fix_event.cpp fix_external.cpp fix_freeze.cpp fix_gravity.cpp fix_heat.cpp fix_indent.cpp fix_langevin.cpp fix_lineforce.cpp fix_minimize.cpp fix_momentum.cpp fix_move.cpp fix_msst.cpp fix_nh.cpp fix_nh_asphere.cpp fix_nh_sphere.cpp fix_nph.cpp fix_nph_asphere.cpp fix_nph_sphere.cpp fix_npt.cpp fix_npt_asphere.cpp fix_npt_sphere.cpp fix_nve.cpp fix_nve_asphere.cpp fix_nve_limit.cpp fix_nve_noforce.cpp fix_nve_sphere.cpp fix_nvt.cpp fix_nvt_asphere.cpp fix_nvt_sllod.cpp fix_nvt_sphere.cpp fix_orient_fcc.cpp fix_peri_neigh.cpp fix_planeforce.cpp fix_pour.cpp fix_press_berendsen.cpp fix_print.cpp fix_qeq_comb.cpp fix_read_restart.cpp fix_recenter.cpp fix_respa.cpp fix_rigid.cpp fix_rigid_nve.cpp fix_rigid_nvt.cpp fix_setforce.cpp fix_shake.cpp fix_shear_history.cpp fix_spring.cpp fix_spring_rg.cpp fix_spring_self.cpp fix_store_force.cpp fix_store_state.cpp fix_temp_berendsen.cpp fix_temp_rescale.cpp fix_thermal_conductivity.cpp fix_tmd.cpp fix_ttm.cpp fix_viscosity.cpp fix_viscous.cpp fix_wall.cpp fix_wall_colloid.cpp fix_wall_gran.cpp fix_wall_harmonic.cpp fix_wall_lj126.cpp fix_wall_lj93.cpp fix_wall_reflect.cpp fix_wall_region.cpp force.cpp group.cpp improper.cpp improper_class2.cpp improper_cvff.cpp improper_harmonic.cpp improper_hybrid.cpp input.cpp integrate.cpp irregular.cpp kspace.cpp lammps.cpp lattice.cpp library.cpp memory.cpp min.cpp min_cg.cpp min_hftn.cpp min_linesearch.cpp min_sd.cpp minimize.cpp modify.cpp neigh_bond.cpp neigh_derive.cpp neigh_full.cpp neigh_gran.cpp neigh_half_bin.cpp neigh_half_multi.cpp neigh_half_nsq.cpp neigh_list.cpp neigh_request.cpp neigh_respa.cpp neigh_stencil.cpp neighbor.cpp output.cpp pack.cpp pair.cpp pair_airebo.cpp pair_born_coul_long.cpp pair_buck.cpp pair_buck_coul_cut.cpp pair_buck_coul_long.cpp pair_colloid.cpp pair_comb.cpp pair_coul_cut.cpp pair_coul_debye.cpp pair_coul_long.cpp pair_dipole_cut.cpp pair_dpd.cpp pair_dpd_tstat.cpp pair_dsmc.cpp pair_eam.cpp pair_eam_alloy.cpp pair_eam_alloy_opt.cpp pair_eam_fs.cpp pair_eam_fs_opt.cpp pair_eam_opt.cpp pair_eim.cpp pair_gayberne.cpp pair_gran_hertz_history.cpp pair_gran_hooke.cpp pair_gran_hooke_history.cpp pair_hybrid.cpp pair_hybrid_overlay.cpp pair_lj96_cut.cpp pair_lj_charmm_coul_charmm.cpp pair_lj_charmm_coul_charmm_implicit.cpp pair_lj_charmm_coul_long.cpp pair_lj_charmm_coul_long_opt.cpp pair_lj_class2.cpp pair_lj_class2_coul_cut.cpp pair_lj_class2_coul_long.cpp pair_lj_cut.cpp pair_lj_cut_coul_cut.cpp pair_lj_cut_coul_debye.cpp pair_lj_cut_coul_long.cpp pair_lj_cut_coul_long_tip4p.cpp pair_lj_cut_opt.cpp pair_lj_expand.cpp pair_lj_gromacs.cpp pair_lj_gromacs_coul_gromacs.cpp pair_lj_smooth.cpp pair_lubricate.cpp pair_morse.cpp pair_morse_opt.cpp pair_peri_lps.cpp pair_peri_pmb.cpp pair_resquared.cpp pair_soft.cpp pair_sw.cpp pair_table.cpp pair_tersoff.cpp pair_tersoff_zbl.cpp pair_yukawa.cpp pair_yukawa_colloid.cpp pppm.cpp pppm_tip4p.cpp prd.cpp random_mars.cpp random_park.cpp read_data.cpp read_restart.cpp region.cpp region_block.cpp region_cone.cpp region_cylinder.cpp region_intersect.cpp region_plane.cpp region_prism.cpp region_sphere.cpp region_union.cpp remap.cpp remap_wrap.cpp replicate.cpp respa.cpp run.cpp set.cpp shell.cpp special.cpp temper.cpp thermo.cpp timer.cpp universe.cpp update.cpp variable.cpp velocity.cpp verlet.cpp write_restart.cpp xdr_compat.cpp -INC = angle.h angle_charmm.h angle_cosine.h angle_cosine_delta.h angle_cosine_squared.h angle_harmonic.h angle_hybrid.h angle_table.h atom.h atom_vec.h atom_vec_angle.h atom_vec_atomic.h atom_vec_bond.h atom_vec_charge.h atom_vec_full.h atom_vec_hybrid.h atom_vec_molecular.h bond.h bond_fene.h bond_fene_expand.h bond_harmonic.h bond_hybrid.h bond_morse.h bond_nonlinear.h bond_quartic.h bond_table.h change_box.h comm.h compute.h compute_angle_local.h compute_bond_local.h compute_centro_atom.h compute_cna_atom.h compute_com.h compute_com_molecule.h compute_coord_atom.h compute_dihedral_local.h compute_displace_atom.h compute_erotate_sphere.h compute_group_group.h compute_gyration.h compute_gyration_molecule.h compute_heat_flux.h compute_improper_local.h compute_ke.h compute_ke_atom.h compute_msd.h compute_msd_molecule.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_pressure.h compute_property_atom.h compute_property_local.h compute_property_molecule.h compute_rdf.h compute_reduce.h compute_reduce_region.h compute_stress_atom.h compute_temp.h compute_temp_com.h compute_temp_deform.h compute_temp_partial.h compute_temp_profile.h compute_temp_ramp.h compute_temp_region.h compute_temp_sphere.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.h dihedral_harmonic.h dihedral_helix.h dihedral_hybrid.h dihedral_multi_harmonic.h dihedral_opls.h displace_atoms.h displace_box.h domain.h dump.h dump_atom.h dump_cfg.h dump_custom.h dump_dcd.h dump_local.h dump_xyz.h error.h ewald.h fft3d.h fft3d_wrap.h finish.h fix.h fix_adapt.h fix_addforce.h fix_ave_atom.h fix_ave_correlate.h fix_ave_histo.h fix_ave_spatial.h fix_ave_time.h fix_aveforce.h fix_bond_break.h fix_bond_create.h fix_bond_swap.h fix_box_relax.h fix_deform.h fix_deposit.h fix_drag.h fix_dt_reset.h fix_efield.h fix_enforce2d.h fix_evaporate.h fix_gravity.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_nh.h fix_nh_sphere.h fix_nph.h fix_nph_sphere.h fix_npt.h fix_npt_sphere.h fix_nve.h fix_nve_limit.h fix_nve_noforce.h fix_nve_sphere.h fix_nvt.h fix_nvt_sllod.h fix_nvt_sphere.h fix_orient_fcc.h fix_planeforce.h fix_press_berendsen.h fix_print.h fix_qeq_comb.h fix_recenter.h fix_respa.h fix_rigid.h fix_rigid_nve.h fix_rigid_nvt.h fix_setforce.h fix_shake.h fix_shear_history.h fix_spring.h fix_spring_rg.h fix_spring_self.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_harmonic.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_reflect.h fix_wall_region.h force.h group.h improper.h improper_cvff.h improper_harmonic.h improper_hybrid.h input.h integrate.h kspace.h lammps.h lattice.h library.h math_extra.h memory.h min.h min_cg.h min_hftn.h min_linesearch.h min_sd.h minimize.h modify.h neigh_list.h neigh_request.h neighbor.h output.h pack.h pair.h pair_airebo.h pair_born_coul_long.h pair_buck.h pair_buck_coul_cut.h pair_buck_coul_long.h pair_comb.h pair_coul_cut.h pair_coul_debye.h pair_coul_long.h pair_dpd.h pair_dpd_tstat.h pair_eam.h pair_eam_alloy.h pair_eam_fs.h pair_eim.h pair_hybrid.h pair_hybrid_overlay.h pair_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.h pair_lj_charmm_coul_long.h pair_lj_cut.h pair_lj_cut_coul_cut.h pair_lj_cut_coul_debye.h pair_lj_cut_coul_long.h pair_lj_cut_coul_long_tip4p.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_smooth.h pair_morse.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_zbl.h pair_yukawa.h pointers.h pppm.h pppm_tip4p.h random_mars.h random_park.h read_data.h read_restart.h region.h region_block.h region_cone.h region_cylinder.h region_intersect.h region_plane.h region_prism.h region_sphere.h region_union.h remap.h remap_wrap.h replicate.h respa.h run.h set.h shell.h special.h style_angle.h style_atom.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_pair.h style_region.h temper.h thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h version.h write_restart.h +INC = angle.h angle_charmm.h angle_class2.h angle_cosine.h angle_cosine_delta.h angle_cosine_squared.h angle_harmonic.h angle_hybrid.h angle_table.h atom.h atom_vec.h atom_vec_angle.h atom_vec_atomic.h atom_vec_bond.h atom_vec_charge.h atom_vec_colloid.h atom_vec_dipole.h atom_vec_ellipsoid.h atom_vec_full.h atom_vec_granular.h atom_vec_hybrid.h atom_vec_molecular.h atom_vec_peri.h bond.h bond_class2.h bond_fene.h bond_fene_expand.h bond_harmonic.h bond_hybrid.h bond_morse.h bond_nonlinear.h bond_quartic.h bond_table.h change_box.h comm.h compute.h compute_angle_local.h compute_bond_local.h compute_centro_atom.h compute_cna_atom.h compute_com.h compute_com_molecule.h compute_coord_atom.h compute_damage_atom.h compute_dihedral_local.h compute_displace_atom.h compute_erotate_asphere.h compute_erotate_sphere.h compute_event_displace.h compute_group_group.h compute_gyration.h compute_gyration_molecule.h compute_heat_flux.h compute_improper_local.h compute_ke.h compute_ke_atom.h compute_msd.h compute_msd_molecule.h compute_pair_local.h compute_pe.h compute_pe_atom.h compute_pressure.h compute_property_atom.h compute_property_local.h compute_property_molecule.h compute_rdf.h compute_reduce.h compute_reduce_region.h compute_stress_atom.h compute_temp.h compute_temp_asphere.h compute_temp_com.h compute_temp_deform.h compute_temp_partial.h compute_temp_profile.h compute_temp_ramp.h compute_temp_region.h compute_temp_sphere.h create_atoms.h create_box.h delete_atoms.h delete_bonds.h dihedral.h dihedral_charmm.h dihedral_class2.h dihedral_harmonic.h dihedral_helix.h dihedral_hybrid.h dihedral_multi_harmonic.h dihedral_opls.h displace_atoms.h displace_box.h domain.h dump.h dump_atom.h dump_cfg.h dump_custom.h dump_dcd.h dump_local.h dump_xtc.h dump_xyz.h error.h ewald.h fft3d.h fft3d_wrap.h finish.h fix.h fix_adapt.h fix_addforce.h fix_ave_atom.h fix_ave_correlate.h fix_ave_histo.h fix_ave_spatial.h fix_ave_time.h fix_aveforce.h fix_bond_break.h fix_bond_create.h fix_bond_swap.h fix_box_relax.h fix_deform.h fix_deposit.h fix_drag.h fix_dt_reset.h fix_efield.h fix_enforce2d.h fix_evaporate.h fix_event.h fix_external.h fix_freeze.h fix_gravity.h fix_heat.h fix_indent.h fix_langevin.h fix_lineforce.h fix_minimize.h fix_momentum.h fix_move.h fix_msst.h fix_nh.h fix_nh_asphere.h fix_nh_sphere.h fix_nph.h fix_nph_asphere.h fix_nph_sphere.h fix_npt.h fix_npt_asphere.h fix_npt_sphere.h fix_nve.h fix_nve_asphere.h fix_nve_limit.h fix_nve_noforce.h fix_nve_sphere.h fix_nvt.h fix_nvt_asphere.h fix_nvt_sllod.h fix_nvt_sphere.h fix_orient_fcc.h fix_peri_neigh.h fix_planeforce.h fix_pour.h fix_press_berendsen.h fix_print.h fix_qeq_comb.h fix_read_restart.h fix_recenter.h fix_respa.h fix_rigid.h fix_rigid_nve.h fix_rigid_nvt.h fix_setforce.h fix_shake.h fix_shear_history.h fix_spring.h fix_spring_rg.h fix_spring_self.h fix_store_force.h fix_store_state.h fix_temp_berendsen.h fix_temp_rescale.h fix_thermal_conductivity.h fix_tmd.h fix_ttm.h fix_viscosity.h fix_viscous.h fix_wall.h fix_wall_colloid.h fix_wall_gran.h fix_wall_harmonic.h fix_wall_lj126.h fix_wall_lj93.h fix_wall_reflect.h fix_wall_region.h force.h group.h improper.h improper_class2.h improper_cvff.h improper_harmonic.h improper_hybrid.h input.h integrate.h irregular.h kspace.h lammps.h lattice.h library.h math_extra.h memory.h min.h min_cg.h min_hftn.h min_linesearch.h min_sd.h minimize.h modify.h neigh_list.h neigh_request.h neighbor.h output.h pack.h pair.h pair_airebo.h pair_born_coul_long.h pair_buck.h pair_buck_coul_cut.h pair_buck_coul_long.h pair_colloid.h pair_comb.h pair_coul_cut.h pair_coul_debye.h pair_coul_long.h pair_dipole_cut.h pair_dpd.h pair_dpd_tstat.h pair_dsmc.h pair_eam.h pair_eam_alloy.h pair_eam_alloy_opt.h pair_eam_fs.h pair_eam_fs_opt.h pair_eam_opt.h pair_eim.h pair_gayberne.h pair_gran_hertz_history.h pair_gran_hooke.h pair_gran_hooke_history.h pair_hybrid.h pair_hybrid_overlay.h pair_lj96_cut.h pair_lj_charmm_coul_charmm.h pair_lj_charmm_coul_charmm_implicit.h pair_lj_charmm_coul_long.h pair_lj_charmm_coul_long_opt.h pair_lj_class2.h pair_lj_class2_coul_cut.h pair_lj_class2_coul_long.h pair_lj_cut.h pair_lj_cut_coul_cut.h pair_lj_cut_coul_debye.h pair_lj_cut_coul_long.h pair_lj_cut_coul_long_tip4p.h pair_lj_cut_opt.h pair_lj_expand.h pair_lj_gromacs.h pair_lj_gromacs_coul_gromacs.h pair_lj_smooth.h pair_lubricate.h pair_morse.h pair_morse_opt.h pair_peri_lps.h pair_peri_pmb.h pair_resquared.h pair_soft.h pair_sw.h pair_table.h pair_tersoff.h pair_tersoff_zbl.h pair_yukawa.h pair_yukawa_colloid.h pointers.h pppm.h pppm_tip4p.h prd.h random_mars.h random_park.h read_data.h read_restart.h region.h region_block.h region_cone.h region_cylinder.h region_intersect.h region_plane.h region_prism.h region_sphere.h region_union.h remap.h remap_wrap.h replicate.h respa.h run.h set.h shell.h special.h style_angle.h style_atom.h style_bond.h style_command.h style_compute.h style_dihedral.h style_dump.h style_fix.h style_improper.h style_integrate.h style_kspace.h style_minimize.h style_pair.h style_region.h temper.h thermo.h timer.h universe.h update.h variable.h velocity.h verlet.h version.h write_restart.h xdr_compat.h OBJ = $(SRC:.cpp=.o)