diff --git a/lib/netcdf/Makefile.lammps b/lib/netcdf/Makefile.lammps new file mode 100644 index 0000000000..3da0bebea4 --- /dev/null +++ b/lib/netcdf/Makefile.lammps @@ -0,0 +1,19 @@ +# Settings that the LAMMPS build will import when a package using the +# netCDF library is installed. This tries to automate configuration +# via the nc-config tool, which is part of the netCDF installation. + +ifeq ($(shell nc-config --has-pnetcdf),no) + +netcdf_SYSINC = $(shell nc-config --cflags) +netcdf_SYSLIB = $(shell nc-config --libs) +netcdf_SYSPATH = + +endif + +ifeq ($(shell nc-config --has-pnetcdf),yes) + +netcdf_SYSINC = -DLMP_HAS_PNETCDF $(nc-config --cflags) +netcdf_SYSLIB = $(shell nc-config --libs) +netcdf_SYSPATH = + +endif diff --git a/src/Makefile b/src/Makefile index 12466e01fa..1e9d3dc86d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -57,7 +57,7 @@ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \ PACKLIB = compress gpu kim kokkos meam mpiio poems python reax voronoi \ user-atc user-awpmd user-colvars user-h5md user-lb user-molfile \ - user-qmmm user-quip user-vtk + user-nc-dump user-qmmm user-quip user-vtk PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/USER-NC-DUMP/Install.sh b/src/USER-NC-DUMP/Install.sh new file mode 100644 index 0000000000..37ebd0a0a5 --- /dev/null +++ b/src/USER-NC-DUMP/Install.sh @@ -0,0 +1,63 @@ +# Install/unInstall package files in LAMMPS +# mode = 0/1/2 for uninstall/install/update + +mode=$1 + +# enforce using portable C locale +LC_ALL=C +export LC_ALL + +# arg1 = file, arg2 = file it depends on + +action () { + if (test $mode = 0) then + rm -f ../$1 + elif (! cmp -s $1 ../$1) then + if (test -z "$2" || test -e ../$2) then + cp $1 .. + if (test $mode = 2) then + echo " updating src/$1" + fi + fi + elif (test -n "$2") then + if (test ! -e ../$2) then + rm -f ../$1 + fi + fi +} + +for file in *.cpp *.h; do + action $file +done + +# edit 2 Makefile.package files to include/exclude package info + +if (test $1 = 1) then + + if (test -e ../Makefile.package) then + sed -i -e 's/[^ \t]*netcdf[^ \t]* //g' ../Makefile.package + sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(netcdf_SYSINC) |' ../Makefile.package + sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(netcdf_SYSLIB) |' ../Makefile.package + sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(netcdf_SYSPATH) |' ../Makefile.package + fi + + if (test -e ../Makefile.package.settings) then + sed -i -e '/^include.*netcdf.*$/d' ../Makefile.package.settings + # multiline form needed for BSD sed on Macs + sed -i -e '4 i \ +include ..\/..\/lib\/netcdf\/Makefile.lammps +' ../Makefile.package.settings + + fi + +elif (test $1 = 0) then + + if (test -e ../Makefile.package) then + sed -i -e 's/[^ \t]*netcdf[^ \t]* //g' ../Makefile.package + fi + + if (test -e ../Makefile.package.settings) then + sed -i -e '/^include.*netcdf.*$/d' ../Makefile.package.settings + fi + +fi diff --git a/src/USER-NC-DUMP/README b/src/USER-NC-DUMP/README new file mode 100644 index 0000000000..fc7c8df121 --- /dev/null +++ b/src/USER-NC-DUMP/README @@ -0,0 +1,37 @@ +USER-NC-DUMP +============ + +This package provides the nc and (optionally) the nc/mpiio dump styles. + +PACKAGE DESCRIPTION +------------------- + +This is a LAMMPS (http://lammps.sandia.gov/) dump style for output into a NetCDF +database. The database format follows the AMBER NetCDF trajectory convention +(http://ambermd.org/netcdf/nctraj.xhtml), but includes extensions to this +convention. These extension are: +* A variable "cell_origin" (of dimension "frame", "cell_spatial") that contains + the bottom left corner of the simulation cell. +* Any number of additional variables corresponding to per atom scalar, vector + or tensor quantities available within LAMMPS. Tensor quantities are written in + Voigt notation. An additional dimension "Voigt" of length 6 is created for + this purpose. +* Possibility to output to an HDF5 database. + +NetCDF files can be directly visualized with the following tools: +* Ovito (http://www.ovito.org/). Ovito supports the AMBER convention and all of + the above extensions. +* VMD (http://www.ks.uiuc.edu/Research/vmd/). +* AtomEye (http://www.libatoms.org/). The libAtoms version of AtomEye contains + a NetCDF reader that is not present in the standard distribution of AtomEye. + +The person who created these files is Lars Pastewka at +Karlsruhe Institute of Technology (lars.pastewka@kit.edu). +Contact him directly if you have questions. + +Lars Pastewka +Institute for Applied Materials (IAM) +Karlsruhe Institute of Technology (KIT) +Kaiserstraße 12, 76131 Karlsruhe +e-mail: lars.pastewka@kit.edu + diff --git a/src/USER-NC-DUMP/dump_nc.cpp b/src/USER-NC-DUMP/dump_nc.cpp new file mode 100644 index 0000000000..5dd483660c --- /dev/null +++ b/src/USER-NC-DUMP/dump_nc.cpp @@ -0,0 +1,1144 @@ +/* ====================================================================== + LAMMPS NetCDF dump style + https://github.com/pastewka/lammps-netcdf + Lars Pastewka, lars.pastewka@kit.edu + + Copyright (2011-2013) Fraunhofer IWM + Copyright (2014) Karlsruhe Institute of Technology + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + ====================================================================== */ + +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include + +#include + +#include "atom.h" +#include "comm.h" +#include "compute.h" +#include "domain.h" +#include "error.h" +#include "fix.h" +#include "group.h" +#include "input.h" +#include "math_const.h" +#include "memory.h" +#include "modify.h" +#include "stdlib.h" +#include "string.h" +#include "update.h" +#include "universe.h" +#include "variable.h" +#include "force.h" + +#include "dump_nc.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +enum{INT,DOUBLE}; // same as in dump_custom.cpp + +const char NC_FRAME_STR[] = "frame"; +const char NC_SPATIAL_STR[] = "spatial"; +const char NC_VOIGT_STR[] = "Voigt"; +const char NC_ATOM_STR[] = "atom"; +const char NC_CELL_SPATIAL_STR[] = "cell_spatial"; +const char NC_CELL_ANGULAR_STR[] = "cell_angular"; +const char NC_LABEL_STR[] = "label"; + +const char NC_TIME_STR[] = "time"; +const char NC_CELL_ORIGIN_STR[] = "cell_origin"; +const char NC_CELL_LENGTHS_STR[] = "cell_lengths"; +const char NC_CELL_ANGLES_STR[] = "cell_angles"; + +const char NC_UNITS_STR[] = "units"; +const char NC_SCALE_FACTOR_STR[] = "scale_factor"; + +const int THIS_IS_A_FIX = -1; +const int THIS_IS_A_COMPUTE = -2; +const int THIS_IS_A_VARIABLE = -3; +const int THIS_IS_A_BIGINT = -4; + +/* ---------------------------------------------------------------------- */ + +#define NCERR(x) ncerr(x, NULL, __LINE__) +#define NCERRX(x, descr) ncerr(x, descr, __LINE__) + +/* ---------------------------------------------------------------------- */ + +DumpNC::DumpNC(LAMMPS *lmp, int narg, char **arg) : + DumpCustom(lmp, narg, arg) +{ + // arrays for data rearrangement + + sort_flag = 1; + sortcol = 0; + binary = 1; + flush_flag = 0; + + if (multiproc) + error->all(FLERR,"Multi-processor writes are not supported."); + if (multifile) + error->all(FLERR,"Multiple files are not supported."); + + perat = new nc_perat_t[nfield]; + + for (int i = 0; i < nfield; i++) { + perat[i].dims = 0; + } + + n_perat = 0; + for (int iarg = 5; iarg < narg; iarg++) { + int i = iarg-5; + int idim = 0; + int ndims = 1; + char mangled[1024]; + bool constant = false; + + strcpy(mangled, arg[iarg]); + + // name mangling + // in the AMBER specification + if (!strcmp(mangled, "x") || !strcmp(mangled, "y") || + !strcmp(mangled, "z")) { + idim = mangled[0] - 'x'; + ndims = 3; + strcpy(mangled, "coordinates"); + } + else if (!strcmp(mangled, "vx") || !strcmp(mangled, "vy") || + !strcmp(mangled, "vz")) { + idim = mangled[1] - 'x'; + ndims = 3; + strcpy(mangled, "velocities"); + } + // extensions to the AMBER specification + else if (!strcmp(mangled, "type")) { + strcpy(mangled, "atom_types"); + } + else if (!strcmp(mangled, "xs") || !strcmp(mangled, "ys") || + !strcmp(mangled, "zs")) { + idim = mangled[0] - 'x'; + ndims = 3; + strcpy(mangled, "scaled_coordinates"); + } + else if (!strcmp(mangled, "xu") || !strcmp(mangled, "yu") || + !strcmp(mangled, "zu")) { + idim = mangled[0] - 'x'; + ndims = 3; + strcpy(mangled, "unwrapped_coordinates"); + } + else if (!strcmp(mangled, "fx") || !strcmp(mangled, "fy") || + !strcmp(mangled, "fz")) { + idim = mangled[1] - 'x'; + ndims = 3; + strcpy(mangled, "forces"); + } + else if (!strcmp(mangled, "mux") || !strcmp(mangled, "muy") || + !strcmp(mangled, "muz")) { + idim = mangled[2] - 'x'; + ndims = 3; + strcpy(mangled, "mu"); + } + else if (!strncmp(mangled, "c_", 2)) { + char *ptr = strchr(mangled, '['); + if (ptr) { + if (mangled[strlen(mangled)-1] != ']') + error->all(FLERR,"Missing ']' in dump command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + ndims = THIS_IS_A_COMPUTE; + } + } + else if (!strncmp(mangled, "f_", 2)) { + char *ptr = strchr(mangled, '['); + if (ptr) { + if (mangled[strlen(mangled)-1] != ']') + error->all(FLERR,"Missing ']' in dump command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + ndims = THIS_IS_A_FIX; + } + } + + // find mangled name + int inc = -1; + for (int j = 0; j < n_perat && inc < 0; j++) { + if (!strcmp(perat[j].name, mangled)) { + inc = j; + } + } + + if (inc < 0) { + // this has not yet been defined + inc = n_perat; + perat[inc].dims = ndims; + if (ndims < 0) ndims = DUMP_NC_MAX_DIMS; + for (int j = 0; j < DUMP_NC_MAX_DIMS; j++) { + perat[inc].field[j] = -1; + } + strcpy(perat[inc].name, mangled); + n_perat++; + } + + perat[inc].constant = constant; + perat[inc].ndumped = 0; + perat[inc].field[idim] = i; + } + + n_perframe = 0; + perframe = NULL; + + n_buffer = 0; + int_buffer = NULL; + double_buffer = NULL; + + double_precision = false; + + framei = 0; +} + +/* ---------------------------------------------------------------------- */ + +DumpNC::~DumpNC() +{ + closefile(); + + delete [] perat; + if (n_perframe > 0) + delete [] perframe; + + if (int_buffer) memory->sfree(int_buffer); + if (double_buffer) memory->sfree(double_buffer); +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::openfile() +{ + // now the computes and fixes have been initialized, so we can query + // for the size of vector quantities + for (int i = 0; i < n_perat; i++) { + if (perat[i].dims == THIS_IS_A_COMPUTE) { + int j = -1; + for (int k = 0; k < DUMP_NC_MAX_DIMS; k++) { + if (perat[i].field[k] >= 0) { + j = field2index[perat[i].field[0]]; + } + } + if (j < 0) + error->all(FLERR,"Internal error."); + if (!compute[j]->peratom_flag) + error->all(FLERR,"compute does not provide per atom data"); + perat[i].dims = compute[j]->size_peratom_cols; + if (perat[i].dims > DUMP_NC_MAX_DIMS) + error->all(FLERR,"perat[i].dims > DUMP_NC_MAX_DIMS"); + } + else if (perat[i].dims == THIS_IS_A_FIX) { + int j = -1; + for (int k = 0; k < DUMP_NC_MAX_DIMS; k++) { + if (perat[i].field[k] >= 0) { + j = field2index[perat[i].field[0]]; + } + } + if (j < 0) + error->all(FLERR,"Internal error."); + if (!fix[j]->peratom_flag) + error->all(FLERR,"fix does not provide per atom data"); + perat[i].dims = fix[j]->size_peratom_cols; + if (perat[i].dims > DUMP_NC_MAX_DIMS) + error->all(FLERR,"perat[i].dims > DUMP_NC_MAX_DIMS"); + } + } + + // get total number of atoms + ntotalgr = group->count(igroup); + + if (filewriter) { + if (append_flag && access(filename, F_OK) != -1) { + // Fixme! Perform checks if dimensions and variables conform with + // data structure standard. + + size_t index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + double d[1]; + + if (singlefile_opened) return; + singlefile_opened = 1; + + NCERRX( nc_open(filename, NC_WRITE, &ncid), filename ); + + // dimensions + NCERRX( nc_inq_dimid(ncid, NC_FRAME_STR, &frame_dim), NC_FRAME_STR ); + NCERRX( nc_inq_dimid(ncid, NC_SPATIAL_STR, &spatial_dim), + NC_SPATIAL_STR ); + NCERRX( nc_inq_dimid(ncid, NC_VOIGT_STR, &Voigt_dim), NC_VOIGT_STR ); + NCERRX( nc_inq_dimid(ncid, NC_ATOM_STR, &atom_dim), NC_ATOM_STR ); + NCERRX( nc_inq_dimid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_dim), + NC_CELL_SPATIAL_STR ); + NCERRX( nc_inq_dimid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_dim), + NC_CELL_ANGULAR_STR ); + NCERRX( nc_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR ); + + // default variables + NCERRX( nc_inq_varid(ncid, NC_SPATIAL_STR, &spatial_var), + NC_SPATIAL_STR ); + NCERRX( nc_inq_varid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_var), + NC_CELL_SPATIAL_STR); + NCERRX( nc_inq_varid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_var), + NC_CELL_ANGULAR_STR); + + NCERRX( nc_inq_varid(ncid, NC_TIME_STR, &time_var), NC_TIME_STR ); + NCERRX( nc_inq_varid(ncid, NC_CELL_ORIGIN_STR, &cell_origin_var), + NC_CELL_ORIGIN_STR ); + NCERRX( nc_inq_varid(ncid, NC_CELL_LENGTHS_STR, &cell_lengths_var), + NC_CELL_LENGTHS_STR); + NCERRX( nc_inq_varid(ncid, NC_CELL_ANGLES_STR, &cell_angles_var), + NC_CELL_ANGLES_STR); + + // variables specified in the input file + for (int i = 0; i < n_perat; i++) { + nc_type xtype; + + // Type mangling + if (vtype[perat[i].field[0]] == INT) { + xtype = NC_INT; + } + else { + if (double_precision) + xtype = NC_DOUBLE; + else + xtype = NC_FLOAT; + } + + NCERRX( nc_inq_varid(ncid, perat[i].name, &perat[i].var), + perat[i].name ); + } + + // perframe variables + for (int i = 0; i < n_perframe; i++) { + NCERRX( nc_inq_varid(ncid, perframe[i].name, &perframe[i].var), + perframe[i].name ); + } + + size_t nframes; + NCERR( nc_inq_dimlen(ncid, frame_dim, &nframes) ); + // framei == -1 means append to file, == -2 means override last frame + // Note that in the input file this translates to 'yes', '-1', etc. + if (framei < 0 || (append_flag && framei == 0)) framei = nframes+framei+1; + if (framei < 1) framei = 1; + } + else { + int dims[NC_MAX_VAR_DIMS]; + size_t index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + double d[1]; + + if (singlefile_opened) return; + singlefile_opened = 1; + + NCERRX( nc_create(filename, NC_64BIT_OFFSET, &ncid), + filename ); + + // dimensions + NCERRX( nc_def_dim(ncid, NC_FRAME_STR, NC_UNLIMITED, &frame_dim), + NC_FRAME_STR ); + NCERRX( nc_def_dim(ncid, NC_SPATIAL_STR, 3, &spatial_dim), + NC_SPATIAL_STR ); + NCERRX( nc_def_dim(ncid, NC_VOIGT_STR, 6, &Voigt_dim), + NC_VOIGT_STR ); + NCERRX( nc_def_dim(ncid, NC_ATOM_STR, ntotalgr, &atom_dim), + NC_ATOM_STR ); + NCERRX( nc_def_dim(ncid, NC_CELL_SPATIAL_STR, 3, &cell_spatial_dim), + NC_CELL_SPATIAL_STR ); + NCERRX( nc_def_dim(ncid, NC_CELL_ANGULAR_STR, 3, &cell_angular_dim), + NC_CELL_ANGULAR_STR ); + NCERRX( nc_def_dim(ncid, NC_LABEL_STR, 10, &label_dim), + NC_LABEL_STR ); + + // default variables + dims[0] = spatial_dim; + NCERRX( nc_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var), + NC_SPATIAL_STR ); + NCERRX( nc_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims, + &cell_spatial_var), NC_CELL_SPATIAL_STR ); + dims[0] = spatial_dim; + dims[1] = label_dim; + NCERRX( nc_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims, + &cell_angular_var), NC_CELL_ANGULAR_STR ); + + dims[0] = frame_dim; + NCERRX( nc_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var), + NC_TIME_STR); + dims[0] = frame_dim; + dims[1] = cell_spatial_dim; + NCERRX( nc_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims, + &cell_origin_var), NC_CELL_ORIGIN_STR ); + NCERRX( nc_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims, + &cell_lengths_var), NC_CELL_LENGTHS_STR ); + dims[0] = frame_dim; + dims[1] = cell_angular_dim; + NCERRX( nc_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims, + &cell_angles_var), NC_CELL_ANGLES_STR ); + + // variables specified in the input file + dims[0] = frame_dim; + dims[1] = atom_dim; + dims[2] = spatial_dim; + + for (int i = 0; i < n_perat; i++) { + nc_type xtype; + + // Type mangling + if (vtype[perat[i].field[0]] == INT) { + xtype = NC_INT; + } + else { + if (double_precision) + xtype = NC_DOUBLE; + else + xtype = NC_FLOAT; + } + + if (perat[i].constant) { + // this quantity will only be written once + if (perat[i].dims == 6) { + // this is a tensor in Voigt notation + dims[2] = Voigt_dim; + NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1, + &perat[i].var), perat[i].name ); + } + else if (perat[i].dims == 3) { + // this is a vector, we need to store x-, y- and z-coordinates + dims[2] = spatial_dim; + NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1, + &perat[i].var), perat[i].name ); + } + else if (perat[i].dims == 1) { + NCERRX( nc_def_var(ncid, perat[i].name, xtype, 1, dims+1, + &perat[i].var), perat[i].name ); + } + else { + char errstr[1024]; + sprintf(errstr, "%i dimensions for '%s'. Not sure how to write " + "this to the NetCDF trajectory file.", perat[i].dims, + perat[i].name); + error->all(FLERR,errstr); + } + } + else { + if (perat[i].dims == 6) { + // this is a tensor in Voigt notation + dims[2] = Voigt_dim; + NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims, + &perat[i].var), perat[i].name ); + } + else if (perat[i].dims == 3) { + // this is a vector, we need to store x-, y- and z-coordinates + dims[2] = spatial_dim; + NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims, + &perat[i].var), perat[i].name ); + } + else if (perat[i].dims == 1) { + NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims, + &perat[i].var), perat[i].name ); + } + else { + char errstr[1024]; + sprintf(errstr, "%i dimensions for '%s'. Not sure how to write " + "this to the NetCDF trajectory file.", perat[i].dims, + perat[i].name); + error->all(FLERR,errstr); + } + } + } + + // perframe variables + for (int i = 0; i < n_perframe; i++) { + if (perframe[i].type == THIS_IS_A_BIGINT) { + NCERRX( nc_def_var(ncid, perframe[i].name, NC_LONG, 1, dims, + &perframe[i].var), perframe[i].name ); + } + else { + NCERRX( nc_def_var(ncid, perframe[i].name, NC_DOUBLE, 1, dims, + &perframe[i].var), perframe[i].name ); + } + } + + // attributes + NCERR( nc_put_att_text(ncid, NC_GLOBAL, "Conventions", + 5, "AMBER") ); + NCERR( nc_put_att_text(ncid, NC_GLOBAL, "ConventionVersion", + 3, "1.0") ); + + NCERR( nc_put_att_text(ncid, NC_GLOBAL, "program", + 6, "LAMMPS") ); + NCERR( nc_put_att_text(ncid, NC_GLOBAL, "programVersion", + strlen(universe->version), universe->version) ); + + // units + if (!strcmp(update->unit_style, "lj")) { + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, + 2, "lj") ); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 2, "lj") ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 2, "lj") ); + } + else if (!strcmp(update->unit_style, "real")) { + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, + 11, "femtosecond") ); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 8, "Angstrom") ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 8, "Angstrom") ); + } + else if (!strcmp(update->unit_style, "metal")) { + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, + 10, "picosecond") ); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 8, "Angstrom") ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 8, "Angstrom") ); + } + else if (!strcmp(update->unit_style, "si")) { + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, + 6, "second") ); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 5, "meter") ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 5, "meter") ); + } + else if (!strcmp(update->unit_style, "cgs")) { + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, + 6, "second") ); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 10, "centimeter") ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 10, "centimeter") ); + } + else if (!strcmp(update->unit_style, "electron")) { + NCERR( nc_put_att_text(ncid, time_var, NC_UNITS_STR, + 11, "femtosecond") ); + NCERR( nc_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 4, "Bohr") ); + NCERR( nc_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 4, "Bohr") ); + } + else { + char errstr[1024]; + sprintf(errstr, "Unsupported unit style '%s'", update->unit_style); + error->all(FLERR,errstr); + } + + NCERR( nc_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, + 6, "degree") ); + + d[0] = update->dt; + NCERR( nc_put_att_double(ncid, time_var, NC_SCALE_FACTOR_STR, + NC_DOUBLE, 1, d) ); + d[0] = 1.0; + NCERR( nc_put_att_double(ncid, cell_origin_var, NC_SCALE_FACTOR_STR, + NC_DOUBLE, 1, d) ); + d[0] = 1.0; + NCERR( nc_put_att_double(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR, + NC_DOUBLE, 1, d) ); + + /* + * Finished with definition + */ + + NCERR( nc_enddef(ncid) ); + + /* + * Write label variables + */ + + NCERR( nc_put_var_text(ncid, spatial_var, "xyz") ); + NCERR( nc_put_var_text(ncid, cell_spatial_var, "abc") ); + index[0] = 0; + index[1] = 0; + count[0] = 1; + count[1] = 5; + NCERR( nc_put_vara_text(ncid, cell_angular_var, index, count, "alpha") ); + index[0] = 1; + count[1] = 4; + NCERR( nc_put_vara_text(ncid, cell_angular_var, index, count, "beta") ); + index[0] = 2; + count[1] = 5; + NCERR( nc_put_vara_text(ncid, cell_angular_var, index, count, "gamma") ); + + framei = 1; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::closefile() +{ + if (filewriter && singlefile_opened) { + NCERR( nc_close(ncid) ); + singlefile_opened = 0; + // append next time DumpNC::openfile is called + append_flag = 1; + // write to next frame upon next open + framei++; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::write() +{ + // open file + + openfile(); + + // need to write per-frame (global) properties here since they may come + // from computes. write_header below is only called from the writing + // processes, but modify->compute[j]->compute_* must be called from all + // processes. + + size_t start[2]; + + start[0] = framei-1; + start[1] = 0; + + for (int i = 0; i < n_perframe; i++) { + + if (perframe[i].type == THIS_IS_A_BIGINT) { + bigint data; + (this->*perframe[i].compute)((void*) &data); + + if (filewriter) +#if defined(LAMMPS_SMALLBIG) || defined(LAMMPS_BIGBIG) + NCERR( nc_put_var1_long(ncid, perframe[i].var, start, &data) ); +#else + NCERR( nc_put_var1_int(ncid, perframe[i].var, start, &data) ); +#endif + } + else { + double data; + int j = perframe[i].index; + int idim = perframe[i].dim; + + if (perframe[i].type == THIS_IS_A_COMPUTE) { + if (idim >= 0) { + modify->compute[j]->compute_vector(); + data = modify->compute[j]->vector[idim]; + } + else + data = modify->compute[j]->compute_scalar(); + } + else if (perframe[i].type == THIS_IS_A_FIX) { + if (idim >= 0) { + data = modify->fix[j]->compute_vector(idim); + } + else + data = modify->fix[j]->compute_scalar(); + } + else if (perframe[i].type == THIS_IS_A_VARIABLE) { + j = input->variable->find(perframe[i].id); + data = input->variable->compute_equal(j); + } + + if (filewriter) + NCERR( nc_put_var1_double(ncid, perframe[i].var, start, &data) ); + } + } + + // call write of superclass + + Dump::write(); + + // close file. this ensures data is flushed and mimized data corruption + + closefile(); +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::write_header(bigint n) +{ + size_t start[2]; + + start[0] = framei-1; + start[1] = 0; + + if (filewriter) { + size_t count[2]; + double time, cell_origin[3], cell_lengths[3], cell_angles[3]; + + time = update->ntimestep; + if (domain->triclinic == 0) { + cell_origin[0] = domain->boxlo[0]; + cell_origin[1] = domain->boxlo[1]; + cell_origin[2] = domain->boxlo[2]; + + cell_lengths[0] = domain->xprd; + cell_lengths[1] = domain->yprd; + cell_lengths[2] = domain->zprd; + + cell_angles[0] = 90; + cell_angles[1] = 90; + cell_angles[2] = 90; + } + else { + double cosalpha, cosbeta, cosgamma; + double *h = domain->h; + + cell_origin[0] = domain->boxlo[0]; + cell_origin[1] = domain->boxlo[1]; + cell_origin[2] = domain->boxlo[2]; + + cell_lengths[0] = domain->xprd; + cell_lengths[1] = sqrt(h[1]*h[1]+h[5]*h[5]); + cell_lengths[2] = sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]); + + cosalpha = (h[5]*h[4]+h[1]*h[3])/ + sqrt((h[1]*h[1]+h[5]*h[5])*(h[2]*h[2]+h[3]*h[3]+h[4]*h[4])); + cosbeta = h[4]/sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]); + cosgamma = h[5]/sqrt(h[1]*h[1]+h[5]*h[5]); + + cell_angles[0] = acos(cosalpha)*180.0/MY_PI; + cell_angles[1] = acos(cosbeta)*180.0/MY_PI; + cell_angles[2] = acos(cosgamma)*180.0/MY_PI; + } + + // Recent AMBER conventions say that nonperiodic boundaries should have + // 'cell_lengths' set to zero. + for (int dim = 0; dim < 3; dim++) { + if (!domain->periodicity[dim]) + cell_lengths[dim] = 0.0; + } + + count[0] = 1; + count[1] = 3; + NCERR( nc_put_var1_double(ncid, time_var, start, &time) ); + NCERR( nc_put_vara_double(ncid, cell_origin_var, start, count, + cell_origin) ); + NCERR( nc_put_vara_double(ncid, cell_lengths_var, start, count, + cell_lengths) ); + NCERR( nc_put_vara_double(ncid, cell_angles_var, start, count, + cell_angles) ); + } + + ndata = n; + blocki = 0; +} + + +/* ---------------------------------------------------------------------- + write data lines to file in a block-by-block style + write head of block (mass & element name) only if has atoms of the type +------------------------------------------------------------------------- */ + +void DumpNC::write_data(int n, double *mybuf) +{ + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + ptrdiff_t stride[NC_MAX_VAR_DIMS]; + + if (!int_buffer) { + n_buffer = n; + int_buffer = (int *) + memory->smalloc(n*sizeof(int), "DumpNC::int_buffer"); + double_buffer = (double *) + memory->smalloc(n*sizeof(double), "DumpNC::double_buffer"); + } + + if (n > n_buffer) { + n_buffer = n; + int_buffer = (int *) + memory->srealloc(int_buffer, n*sizeof(int), "DumpNC::int_buffer"); + double_buffer = (double *) + memory->srealloc(double_buffer, n*sizeof(double), + "DumpNC::double_buffer"); + } + + start[0] = framei-1; + start[1] = blocki; + start[2] = 0; + + count[0] = 1; + count[1] = n; + count[2] = 1; + + stride[0] = 1; + stride[1] = 1; + stride[2] = 3; + + for (int i = 0; i < n_perat; i++) { + int iaux = perat[i].field[0]; + + if (vtype[iaux] == INT) { + // integers + if (perat[i].dims > 1) { + + for (int idim = 0; idim < perat[i].dims; idim++) { + iaux = perat[i].field[idim]; + + if (iaux >= 0) { + for (int j = 0; j < n; j++, iaux+=size_one) { + int_buffer[j] = mybuf[iaux]; + } + + start[2] = idim; + + if (perat[i].constant) { + if (perat[i].ndumped < ntotalgr) { + NCERR( nc_put_vars_int(ncid, perat[i].var, + start+1, count+1, stride+1, + int_buffer) ); + perat[i].ndumped += n; + } + } + else + NCERR( nc_put_vars_int(ncid, perat[i].var, start, count, stride, + int_buffer) ); + } + } + } + else { + for (int j = 0; j < n; j++, iaux+=size_one) { + int_buffer[j] = mybuf[iaux]; + } + + if (perat[i].constant) { + if (perat[i].ndumped < ntotalgr) { + NCERR( nc_put_vara_int(ncid, perat[i].var, start+1, count+1, + int_buffer) ); + perat[i].ndumped += n; + } + } + else + NCERR( nc_put_vara_int(ncid, perat[i].var, start, count, + int_buffer) ); + } + } + else { + // doubles + if (perat[i].dims > 1) { + + for (int idim = 0; idim < perat[i].dims; idim++) { + iaux = perat[i].field[idim]; + + if (iaux >= 0) { + for (int j = 0; j < n; j++, iaux+=size_one) { + double_buffer[j] = mybuf[iaux]; + } + + start[2] = idim; + + if (perat[i].constant) { + if (perat[i].ndumped < ntotalgr) { + NCERR( nc_put_vars_double(ncid, perat[i].var, + start+1, count+1, stride+1, + double_buffer) ); + perat[i].ndumped += n; + } + } + else + NCERR( nc_put_vars_double(ncid, perat[i].var, start, count, + stride, double_buffer) ); + } + } + } + else { + for (int j = 0; j < n; j++, iaux+=size_one) { + double_buffer[j] = mybuf[iaux]; + } + + if (perat[i].constant) { + if (perat[i].ndumped < ntotalgr) { + NCERR( nc_put_vara_double(ncid, perat[i].var, start+1, count+1, + double_buffer) ); + perat[i].ndumped += n; + } + } + else + NCERR( nc_put_vara_double(ncid, perat[i].var, start, count, + double_buffer) ); + } + } + } + + blocki += n; +} + +/* ---------------------------------------------------------------------- */ + +int DumpNC::modify_param(int narg, char **arg) +{ + int iarg = 0; + if (strcmp(arg[iarg],"double") == 0) { + iarg++; + if (iarg >= narg) + error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword."); + if (strcmp(arg[iarg],"yes") == 0) { + double_precision = true; + } + else if (strcmp(arg[iarg],"no") == 0) { + double_precision = false; + } + else error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword."); + iarg++; + return 2; + } + else if (strcmp(arg[iarg],"at") == 0) { + iarg++; + framei = force->inumeric(FLERR,arg[iarg]); + if (framei < 0) framei--; + iarg++; + return 2; + } + else if (strcmp(arg[iarg],"global") == 0) { + // "perframe" quantities, i.e. not per-atom stuff + + iarg++; + + n_perframe = narg-iarg; + perframe = new nc_perframe_t[n_perframe]; + + for (int i = 0; iarg < narg; iarg++, i++) { + int n; + char *suffix; + + if (!strcmp(arg[iarg],"step")) { + perframe[i].type = THIS_IS_A_BIGINT; + perframe[i].compute = &DumpNC::compute_step; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strcmp(arg[iarg],"elapsed")) { + perframe[i].type = THIS_IS_A_BIGINT; + perframe[i].compute = &DumpNC::compute_elapsed; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strcmp(arg[iarg],"elaplong")) { + perframe[i].type = THIS_IS_A_BIGINT; + perframe[i].compute = &DumpNC::compute_elapsed_long; + strcpy(perframe[i].name, arg[iarg]); + } + else { + + n = strlen(arg[iarg]); + + if (n > 2) { + suffix = new char[n-1]; + strcpy(suffix, arg[iarg]+2); + } + else { + char errstr[1024]; + sprintf(errstr, "perframe quantity '%s' must thermo quantity or " + "compute, fix or variable", arg[iarg]); + error->all(FLERR,errstr); + } + + if (!strncmp(arg[iarg], "c_", 2)) { + int idim = -1; + char *ptr = strchr(suffix, '['); + + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all(FLERR,"Missing ']' in dump modify command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + } + + n = modify->find_compute(suffix); + if (n < 0) + error->all(FLERR,"Could not find dump modify compute ID"); + if (modify->compute[n]->peratom_flag != 0) + error->all(FLERR,"Dump modify compute ID computes per-atom info"); + if (idim >= 0 && modify->compute[n]->vector_flag == 0) + error->all(FLERR,"Dump modify compute ID does not compute vector"); + if (idim < 0 && modify->compute[n]->scalar_flag == 0) + error->all(FLERR,"Dump modify compute ID does not compute scalar"); + + perframe[i].type = THIS_IS_A_COMPUTE; + perframe[i].dim = idim; + perframe[i].index = n; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strncmp(arg[iarg], "f_", 2)) { + int idim = -1; + char *ptr = strchr(suffix, '['); + + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all(FLERR,"Missing ']' in dump modify command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + } + + n = modify->find_fix(suffix); + if (n < 0) + error->all(FLERR,"Could not find dump modify fix ID"); + if (modify->fix[n]->peratom_flag != 0) + error->all(FLERR,"Dump modify fix ID computes per-atom info"); + if (idim >= 0 && modify->fix[n]->vector_flag == 0) + error->all(FLERR,"Dump modify fix ID does not compute vector"); + if (idim < 0 && modify->fix[n]->scalar_flag == 0) + error->all(FLERR,"Dump modify fix ID does not compute vector"); + + perframe[i].type = THIS_IS_A_FIX; + perframe[i].dim = idim; + perframe[i].index = n; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strncmp(arg[iarg], "v_", 2)) { + n = input->variable->find(suffix); + if (n < 0) + error->all(FLERR,"Could not find dump modify variable ID"); + if (!input->variable->equalstyle(n)) + error->all(FLERR,"Dump modify variable must be of style equal"); + + perframe[i].type = THIS_IS_A_VARIABLE; + perframe[i].dim = 1; + perframe[i].index = n; + strcpy(perframe[i].name, arg[iarg]); + strcpy(perframe[i].id, suffix); + } + else { + char errstr[1024]; + sprintf(errstr, "perframe quantity '%s' must be compute, fix or " + "variable", arg[iarg]); + error->all(FLERR,errstr); + } + + delete [] suffix; + + } + } + + return narg; + } else return 0; +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::write_prmtop() +{ + char fn[1024]; + char tmp[81]; + FILE *f; + + strcpy(fn, filename); + strcat(fn, ".prmtop"); + + f = fopen(fn, "w"); + fprintf(f, "%%VERSION LAMMPS\n"); + fprintf(f, "%%FLAG TITLE\n"); + fprintf(f, "%%FORMAT(20a4)\n"); + memset(tmp, ' ', 76); + tmp[76] = '\0'; + fprintf(f, "NASN%s\n", tmp); + + fprintf(f, "%%FLAG POINTERS\n"); + fprintf(f, "%%FORMAT(10I8)\n"); +#if defined(LAMMPS_SMALLBIG) || defined(LAMMPS_BIGBIG) + fprintf(f, "%8li", ntotalgr); +#else + fprintf(f, "%8i", ntotalgr); +#endif + for (int i = 0; i < 11; i++) + fprintf(f, "%8i", 0); + fprintf(f, "\n"); + for (int i = 0; i < 12; i++) + fprintf(f, "%8i", 0); + fprintf(f, "\n"); + for (int i = 0; i < 6; i++) + fprintf(f, "%8i", 0); + fprintf(f, "\n"); + + fprintf(f, "%%FLAG ATOM_NAME\n"); + fprintf(f, "%%FORMAT(20a4)\n"); + for (int i = 0; i < ntotalgr; i++) { + fprintf(f, "%4s", "He"); + if ((i+1) % 20 == 0) + fprintf(f, "\n"); + } + + fprintf(f, "%%FLAG CHARGE\n"); + fprintf(f, "%%FORMAT(5E16.5)\n"); + for (int i = 0; i < ntotalgr; i++) { + fprintf(f, "%16.5e", 0.0); + if ((i+1) % 5 == 0) + fprintf(f, "\n"); + } + + fprintf(f, "%%FLAG MASS\n"); + fprintf(f, "%%FORMAT(5E16.5)\n"); + for (int i = 0; i < ntotalgr; i++) { + fprintf(f, "%16.5e", 1.0); + if ((i+1) % 5 == 0) + fprintf(f, "\n"); + } + fclose(f); +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::ncerr(int err, const char *descr, int line) +{ + if (err != NC_NOERR) { + char errstr[1024]; + if (descr) { + sprintf(errstr, "NetCDF failed with error '%s' (while accessing '%s') " + " in line %i of %s.", nc_strerror(err), descr, line, __FILE__); + } + else { + sprintf(errstr, "NetCDF failed with error '%s' in line %i of %s.", + nc_strerror(err), line, __FILE__); + } + error->one(FLERR,errstr); + } +} + +/* ---------------------------------------------------------------------- + one method for every keyword thermo can output + called by compute() or evaluate_keyword() + compute will have already been called + set ivalue/dvalue/bivalue if value is int/double/bigint + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +void DumpNC::compute_step(void *r) +{ + *((bigint *) r) = update->ntimestep; +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::compute_elapsed(void *r) +{ + *((bigint *) r) = update->ntimestep - update->firststep; +} + +/* ---------------------------------------------------------------------- */ + +void DumpNC::compute_elapsed_long(void *r) +{ + *((bigint *) r) = update->ntimestep - update->beginstep; +} diff --git a/src/USER-NC-DUMP/dump_nc.h b/src/USER-NC-DUMP/dump_nc.h new file mode 100644 index 0000000000..75236c15a1 --- /dev/null +++ b/src/USER-NC-DUMP/dump_nc.h @@ -0,0 +1,141 @@ +/* ====================================================================== + LAMMPS NetCDF dump style + https://github.com/pastewka/lammps-netcdf + Lars Pastewka, lars.pastewka@kit.edu + + Copyright (2011-2013) Fraunhofer IWM + Copyright (2014) Karlsruhe Institute of Technology + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + ====================================================================== */ + +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef DUMP_CLASS + +DumpStyle(nc,DumpNC) + +#else + +#ifndef LMP_DUMP_NC_H +#define LMP_DUMP_NC_H + +#include "dump_custom.h" + +namespace LAMMPS_NS { + +const int NC_FIELD_NAME_MAX = 100; +const int DUMP_NC_MAX_DIMS = 100; + +class DumpNC : public DumpCustom { + public: + DumpNC(class LAMMPS *, int, char **); + virtual ~DumpNC(); + virtual void write(); + + private: + // per-atoms quantities (positions, velocities, etc.) + struct nc_perat_t { + int dims; // number of dimensions + int field[DUMP_NC_MAX_DIMS]; // field indices corresponding to the dim. + char name[NC_FIELD_NAME_MAX]; // field name + int var; // NetCDF variable + + bool constant; // is this property per file (not per frame) + int ndumped; // number of enties written for this prop. + }; + + typedef void (DumpNC::*funcptr_t)(void *); + + // per-frame quantities (variables, fixes or computes) + struct nc_perframe_t { + char name[NC_FIELD_NAME_MAX]; // field name + int var; // NetCDF variable + int type; // variable, fix, compute or callback + int index; // index in fix/compute list + funcptr_t compute; // compute function + int dim; // dimension + char id[NC_FIELD_NAME_MAX]; // variable id + + bigint bigint_data; // actual data + double double_data; // actual data + }; + + int framei; // current frame index + int blocki; // current block index + int ndata; // number of data blocks to expect + + bigint ntotalgr; // # of atoms + + int n_perat; // # of netcdf per-atom properties + nc_perat_t *perat; // per-atom properties + + int n_perframe; // # of global netcdf (not per-atom) fix props + nc_perframe_t *perframe; // global properties + + bool double_precision; // write everything as double precision + + bigint n_buffer; // size of buffer + int *int_buffer; // buffer for passing data to netcdf + double *double_buffer; // buffer for passing data to netcdf + + int ncid; + + int frame_dim; + int spatial_dim; + int Voigt_dim; + int atom_dim; + int cell_spatial_dim; + int cell_angular_dim; + int label_dim; + + int spatial_var; + int cell_spatial_var; + int cell_angular_var; + + int time_var; + int cell_origin_var; + int cell_lengths_var; + int cell_angles_var; + + virtual void openfile(); + void closefile(); + virtual void write_header(bigint); + virtual void write_data(int, double *); + void write_prmtop(); + + virtual int modify_param(int, char **); + + void ncerr(int, const char *, int); + + void compute_step(void *); + void compute_elapsed(void *); + void compute_elapsed_long(void *); +}; + +} + +#endif +#endif diff --git a/src/USER-NC-DUMP/dump_nc_mpiio.cpp b/src/USER-NC-DUMP/dump_nc_mpiio.cpp new file mode 100644 index 0000000000..f8960e6901 --- /dev/null +++ b/src/USER-NC-DUMP/dump_nc_mpiio.cpp @@ -0,0 +1,1077 @@ +/* ====================================================================== + LAMMPS NetCDF dump style + https://github.com/pastewka/lammps-netcdf + Lars Pastewka, lars.pastewka@kit.edu + + Copyright (2011-2013) Fraunhofer IWM + Copyright (2014) Karlsruhe Institute of Technology + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + ====================================================================== */ + +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +#if defined(LMP_HAS_PNETCDF) + +#include + +#include + +#include "atom.h" +#include "comm.h" +#include "compute.h" +#include "domain.h" +#include "error.h" +#include "fix.h" +#include "group.h" +#include "input.h" +#include "math_const.h" +#include "memory.h" +#include "modify.h" +#include "stdlib.h" +#include "string.h" +#include "update.h" +#include "universe.h" +#include "variable.h" +#include "force.h" + +#include "dump_nc_mpiio.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +enum{INT,DOUBLE}; // same as in dump_custom.cpp + +const char NC_FRAME_STR[] = "frame"; +const char NC_SPATIAL_STR[] = "spatial"; +const char NC_VOIGT_STR[] = "Voigt"; +const char NC_ATOM_STR[] = "atom"; +const char NC_CELL_SPATIAL_STR[] = "cell_spatial"; +const char NC_CELL_ANGULAR_STR[] = "cell_angular"; +const char NC_LABEL_STR[] = "label"; + +const char NC_TIME_STR[] = "time"; +const char NC_CELL_ORIGIN_STR[] = "cell_origin"; +const char NC_CELL_LENGTHS_STR[] = "cell_lengths"; +const char NC_CELL_ANGLES_STR[] = "cell_angles"; + +const char NC_UNITS_STR[] = "units"; +const char NC_SCALE_FACTOR_STR[] = "scale_factor"; + +const int THIS_IS_A_FIX = -1; +const int THIS_IS_A_COMPUTE = -2; +const int THIS_IS_A_VARIABLE = -3; +const int THIS_IS_A_BIGINT = -4; + +/* ---------------------------------------------------------------------- */ + +#define NCERR(x) ncerr(x, NULL, __LINE__) +#define NCERRX(x, descr) ncerr(x, descr, __LINE__) + +/* ---------------------------------------------------------------------- */ + +DumpNCMPIIO::DumpNCMPIIO(LAMMPS *lmp, int narg, char **arg) : + DumpCustom(lmp, narg, arg) +{ + // arrays for data rearrangement + + sort_flag = 1; + sortcol = 0; + binary = 1; + flush_flag = 0; + + if (multiproc) + error->all(FLERR,"Multi-processor writes are not supported."); + if (multifile) + error->all(FLERR,"Multiple files are not supported."); + + perat = new nc_perat_t[nfield]; + + for (int i = 0; i < nfield; i++) { + perat[i].dims = 0; + } + + n_perat = 0; + for (int iarg = 5; iarg < narg; iarg++) { + int i = iarg-5; + int idim = 0; + int ndims = 1; + char mangled[1024]; + + strcpy(mangled, arg[iarg]); + + // name mangling + // in the AMBER specification + if (!strcmp(mangled, "x") || !strcmp(mangled, "y") || + !strcmp(mangled, "z")) { + idim = mangled[0] - 'x'; + ndims = 3; + strcpy(mangled, "coordinates"); + } + else if (!strcmp(mangled, "vx") || !strcmp(mangled, "vy") || + !strcmp(mangled, "vz")) { + idim = mangled[1] - 'x'; + ndims = 3; + strcpy(mangled, "velocities"); + } + else if (!strcmp(mangled, "xs") || !strcmp(mangled, "ys") || + !strcmp(mangled, "zs")) { + idim = mangled[0] - 'x'; + ndims = 3; + strcpy(mangled, "scaled_coordinates"); + } + else if (!strcmp(mangled, "xu") || !strcmp(mangled, "yu") || + !strcmp(mangled, "zu")) { + idim = mangled[0] - 'x'; + ndims = 3; + strcpy(mangled, "unwrapped_coordinates"); + } + else if (!strcmp(mangled, "fx") || !strcmp(mangled, "fy") || + !strcmp(mangled, "fz")) { + idim = mangled[1] - 'x'; + ndims = 3; + strcpy(mangled, "forces"); + } + else if (!strcmp(mangled, "mux") || !strcmp(mangled, "muy") || + !strcmp(mangled, "muz")) { + idim = mangled[2] - 'x'; + ndims = 3; + strcpy(mangled, "mu"); + } + else if (!strncmp(mangled, "c_", 2)) { + char *ptr = strchr(mangled, '['); + if (ptr) { + if (mangled[strlen(mangled)-1] != ']') + error->all(FLERR,"Missing ']' in dump command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + ndims = THIS_IS_A_COMPUTE; + } + } + else if (!strncmp(mangled, "f_", 2)) { + char *ptr = strchr(mangled, '['); + if (ptr) { + if (mangled[strlen(mangled)-1] != ']') + error->all(FLERR,"Missing ']' in dump command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + ndims = THIS_IS_A_FIX; + } + } + + // find mangled name + int inc = -1; + for (int j = 0; j < n_perat && inc < 0; j++) { + if (!strcmp(perat[j].name, mangled)) { + inc = j; + } + } + + if (inc < 0) { + // this has not yet been defined + inc = n_perat; + perat[inc].dims = ndims; + if (ndims < 0) ndims = DUMP_NC_MPIIO_MAX_DIMS; + for (int j = 0; j < DUMP_NC_MPIIO_MAX_DIMS; j++) { + perat[inc].field[j] = -1; + } + strcpy(perat[inc].name, mangled); + n_perat++; + } + + perat[inc].field[idim] = i; + } + + n_perframe = 0; + perframe = NULL; + + n_buffer = 0; + int_buffer = NULL; + double_buffer = NULL; + + double_precision = false; + + framei = 0; +} + +/* ---------------------------------------------------------------------- */ + +DumpNCMPIIO::~DumpNCMPIIO() +{ + closefile(); + + delete [] perat; + if (n_perframe > 0) + delete [] perframe; + + if (int_buffer) memory->sfree(int_buffer); + if (double_buffer) memory->sfree(double_buffer); +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::openfile() +{ + // now the computes and fixes have been initialized, so we can query + // for the size of vector quantities + for (int i = 0; i < n_perat; i++) { + if (perat[i].dims == THIS_IS_A_COMPUTE) { + int j = -1; + for (int k = 0; k < DUMP_NC_MPIIO_MAX_DIMS; k++) { + if (perat[i].field[k] >= 0) { + j = field2index[perat[i].field[0]]; + } + } + if (j < 0) + error->all(FLERR,"Internal error."); + if (!compute[j]->peratom_flag) + error->all(FLERR,"compute does not provide per atom data"); + perat[i].dims = compute[j]->size_peratom_cols; + if (perat[i].dims > DUMP_NC_MPIIO_MAX_DIMS) + error->all(FLERR,"perat[i].dims > DUMP_NC_MPIIO_MAX_DIMS"); + } + else if (perat[i].dims == THIS_IS_A_FIX) { + int j = -1; + for (int k = 0; k < DUMP_NC_MPIIO_MAX_DIMS; k++) { + if (perat[i].field[k] >= 0) { + j = field2index[perat[i].field[0]]; + } + } + if (j < 0) + error->all(FLERR,"Internal error."); + if (!fix[j]->peratom_flag) + error->all(FLERR,"fix does not provide per atom data"); + perat[i].dims = fix[j]->size_peratom_cols; + if (perat[i].dims > DUMP_NC_MPIIO_MAX_DIMS) + error->all(FLERR,"perat[i].dims > DUMP_NC_MPIIO_MAX_DIMS"); + } + } + + // get total number of atoms + ntotalgr = group->count(igroup); + + if (append_flag && access(filename, F_OK) != -1) { + // Fixme! Perform checks if dimensions and variables conform with + // data structure standard. + + MPI_Offset index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + double d[1]; + + if (singlefile_opened) return; + singlefile_opened = 1; + + NCERRX( ncmpi_open(MPI_COMM_WORLD, filename, NC_WRITE, MPI_INFO_NULL, + &ncid), filename ); + + // dimensions + NCERRX( ncmpi_inq_dimid(ncid, NC_FRAME_STR, &frame_dim), NC_FRAME_STR ); + NCERRX( ncmpi_inq_dimid(ncid, NC_SPATIAL_STR, &spatial_dim), + NC_SPATIAL_STR ); + NCERRX( ncmpi_inq_dimid(ncid, NC_VOIGT_STR, &Voigt_dim), NC_VOIGT_STR ); + NCERRX( ncmpi_inq_dimid(ncid, NC_ATOM_STR, &atom_dim), NC_ATOM_STR ); + NCERRX( ncmpi_inq_dimid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_dim), + NC_CELL_SPATIAL_STR ); + NCERRX( ncmpi_inq_dimid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_dim), + NC_CELL_ANGULAR_STR ); + NCERRX( ncmpi_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR ); + + // default variables + NCERRX( ncmpi_inq_varid(ncid, NC_SPATIAL_STR, &spatial_var), + NC_SPATIAL_STR ); + NCERRX( ncmpi_inq_varid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_var), + NC_CELL_SPATIAL_STR); + NCERRX( ncmpi_inq_varid(ncid, NC_CELL_ANGULAR_STR, &cell_angular_var), + NC_CELL_ANGULAR_STR); + + NCERRX( ncmpi_inq_varid(ncid, NC_TIME_STR, &time_var), NC_TIME_STR ); + NCERRX( ncmpi_inq_varid(ncid, NC_CELL_ORIGIN_STR, &cell_origin_var), + NC_CELL_ORIGIN_STR ); + NCERRX( ncmpi_inq_varid(ncid, NC_CELL_LENGTHS_STR, &cell_lengths_var), + NC_CELL_LENGTHS_STR); + NCERRX( ncmpi_inq_varid(ncid, NC_CELL_ANGLES_STR, &cell_angles_var), + NC_CELL_ANGLES_STR); + + // variables specified in the input file + for (int i = 0; i < n_perat; i++) { + nc_type xtype; + + // Type mangling + if (vtype[perat[i].field[0]] == INT) { + xtype = NC_INT; + } + else { + if (double_precision) + xtype = NC_DOUBLE; + else + xtype = NC_FLOAT; + } + + NCERRX( ncmpi_inq_varid(ncid, perat[i].name, &perat[i].var), + perat[i].name ); + } + + // perframe variables + for (int i = 0; i < n_perframe; i++) { + NCERRX( ncmpi_inq_varid(ncid, perframe[i].name, &perframe[i].var), + perframe[i].name ); + } + + MPI_Offset nframes; + NCERR( ncmpi_inq_dimlen(ncid, frame_dim, &nframes) ); + // framei == -1 means append to file, == -2 means override last frame + // Note that in the input file this translates to 'yes', '-1', etc. + if (framei < 0 || (append_flag && framei == 0)) framei = nframes+framei+1; + if (framei < 1) framei = 1; + } + else { + int dims[NC_MAX_VAR_DIMS]; + MPI_Offset index[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + double d[1]; + + if (singlefile_opened) return; + singlefile_opened = 1; + + NCERRX( ncmpi_create(MPI_COMM_WORLD, filename, NC_64BIT_OFFSET, + MPI_INFO_NULL, &ncid), filename ); + + // dimensions + NCERRX( ncmpi_def_dim(ncid, NC_FRAME_STR, NC_UNLIMITED, &frame_dim), + NC_FRAME_STR ); + NCERRX( ncmpi_def_dim(ncid, NC_SPATIAL_STR, 3, &spatial_dim), + NC_SPATIAL_STR ); + NCERRX( ncmpi_def_dim(ncid, NC_VOIGT_STR, 6, &Voigt_dim), + NC_VOIGT_STR ); + NCERRX( ncmpi_def_dim(ncid, NC_ATOM_STR, ntotalgr, &atom_dim), + NC_ATOM_STR ); + NCERRX( ncmpi_def_dim(ncid, NC_CELL_SPATIAL_STR, 3, &cell_spatial_dim), + NC_CELL_SPATIAL_STR ); + NCERRX( ncmpi_def_dim(ncid, NC_CELL_ANGULAR_STR, 3, &cell_angular_dim), + NC_CELL_ANGULAR_STR ); + NCERRX( ncmpi_def_dim(ncid, NC_LABEL_STR, 10, &label_dim), + NC_LABEL_STR ); + + // default variables + dims[0] = spatial_dim; + NCERRX( ncmpi_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var), + NC_SPATIAL_STR ); + NCERRX( ncmpi_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims, + &cell_spatial_var), NC_CELL_SPATIAL_STR ); + dims[0] = spatial_dim; + dims[1] = label_dim; + NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims, + &cell_angular_var), NC_CELL_ANGULAR_STR ); + + dims[0] = frame_dim; + NCERRX( ncmpi_def_var(ncid, NC_TIME_STR, NC_DOUBLE, 1, dims, &time_var), + NC_TIME_STR); + dims[0] = frame_dim; + dims[1] = cell_spatial_dim; + NCERRX( ncmpi_def_var(ncid, NC_CELL_ORIGIN_STR, NC_DOUBLE, 2, dims, + &cell_origin_var), NC_CELL_ORIGIN_STR ); + NCERRX( ncmpi_def_var(ncid, NC_CELL_LENGTHS_STR, NC_DOUBLE, 2, dims, + &cell_lengths_var), NC_CELL_LENGTHS_STR ); + dims[0] = frame_dim; + dims[1] = cell_angular_dim; + NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGLES_STR, NC_DOUBLE, 2, dims, + &cell_angles_var), NC_CELL_ANGLES_STR ); + + // variables specified in the input file + dims[0] = frame_dim; + dims[1] = atom_dim; + dims[2] = spatial_dim; + + for (int i = 0; i < n_perat; i++) { + nc_type xtype; + + // Type mangling + if (vtype[perat[i].field[0]] == INT) { + xtype = NC_INT; + } + else { + if (double_precision) + xtype = NC_DOUBLE; + else + xtype = NC_FLOAT; + } + + if (perat[i].dims == 6) { + // this is a tensor in Voigt notation + dims[2] = Voigt_dim; + NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims, + &perat[i].var), perat[i].name ); + } + else if (perat[i].dims == 3) { + // this is a vector, we need to store x-, y- and z-coordinates + dims[2] = spatial_dim; + NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims, + &perat[i].var), perat[i].name ); + } + else if (perat[i].dims == 1) { + NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 2, dims, + &perat[i].var), perat[i].name ); + } + else { + char errstr[1024]; + sprintf(errstr, "%i dimensions for '%s'. Not sure how to write " + "this to the NetCDF trajectory file.", perat[i].dims, + perat[i].name); + error->all(FLERR,errstr); + } + } + + // perframe variables + for (int i = 0; i < n_perframe; i++) { + if (perframe[i].type == THIS_IS_A_BIGINT) { + NCERRX( ncmpi_def_var(ncid, perframe[i].name, NC_INT, 1, dims, + &perframe[i].var), perframe[i].name ); + } + else { + NCERRX( ncmpi_def_var(ncid, perframe[i].name, NC_DOUBLE, 1, dims, + &perframe[i].var), perframe[i].name ); + } + } + + // attributes + NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "Conventions", + 5, "AMBER") ); + NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "ConventionVersion", + 3, "1.0") ); + + NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "program", + 6, "LAMMPS") ); + NCERR( ncmpi_put_att_text(ncid, NC_GLOBAL, "programVersion", + strlen(universe->version), universe->version) ); + + // units + if (!strcmp(update->unit_style, "lj")) { + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, + 2, "lj") ); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 2, "lj") ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 2, "lj") ); + } + else if (!strcmp(update->unit_style, "real")) { + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, + 11, "femtosecond") ); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 8, "Angstrom") ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 8, "Angstrom") ); + } + else if (!strcmp(update->unit_style, "metal")) { + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, + 10, "picosecond") ); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 8, "Angstrom") ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 8, "Angstrom") ); + } + else if (!strcmp(update->unit_style, "si")) { + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, + 6, "second") ); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 5, "meter") ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 5, "meter") ); + } + else if (!strcmp(update->unit_style, "cgs")) { + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, + 6, "second") ); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 10, "centimeter") ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 10, "centimeter") ); + } + else if (!strcmp(update->unit_style, "electron")) { + NCERR( ncmpi_put_att_text(ncid, time_var, NC_UNITS_STR, + 11, "femtosecond") ); + NCERR( ncmpi_put_att_text(ncid, cell_origin_var, NC_UNITS_STR, + 4, "Bohr") ); + NCERR( ncmpi_put_att_text(ncid, cell_lengths_var, NC_UNITS_STR, + 4, "Bohr") ); + } + else { + char errstr[1024]; + sprintf(errstr, "Unsupported unit style '%s'", update->unit_style); + error->all(FLERR,errstr); + } + + NCERR( ncmpi_put_att_text(ncid, cell_angles_var, NC_UNITS_STR, + 6, "degree") ); + + d[0] = update->dt; + NCERR( ncmpi_put_att_double(ncid, time_var, NC_SCALE_FACTOR_STR, + NC_DOUBLE, 1, d) ); + d[0] = 1.0; + NCERR( ncmpi_put_att_double(ncid, cell_origin_var, NC_SCALE_FACTOR_STR, + NC_DOUBLE, 1, d) ); + d[0] = 1.0; + NCERR( ncmpi_put_att_double(ncid, cell_lengths_var, NC_SCALE_FACTOR_STR, + NC_DOUBLE, 1, d) ); + + /* + * Finished with definition + */ + + NCERR( ncmpi_enddef(ncid) ); + + /* + * Write label variables + */ + + NCERR( ncmpi_begin_indep_data(ncid) ); + + if (filewriter) { + NCERR( ncmpi_put_var_text(ncid, spatial_var, "xyz") ); + NCERR( ncmpi_put_var_text(ncid, cell_spatial_var, "abc") ); + index[0] = 0; + index[1] = 0; + count[0] = 1; + count[1] = 5; + NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, + "alpha") ); + index[0] = 1; + count[1] = 4; + NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, + "beta") ); + index[0] = 2; + count[1] = 5; + NCERR( ncmpi_put_vara_text(ncid, cell_angular_var, index, count, + "gamma") ); + } + + NCERR( ncmpi_end_indep_data(ncid) ); + + framei = 1; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::closefile() +{ + if (singlefile_opened) { + NCERR( ncmpi_close(ncid) ); + singlefile_opened = 0; + // append next time DumpNCMPIIO::openfile is called + append_flag = 1; + // write to next frame upon next open + framei++; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::write() +{ + // open file + + openfile(); + + // need to write per-frame (global) properties here since they may come + // from computes. write_header below is only called from the writing + // processes, but modify->compute[j]->compute_* must be called from all + // processes. + + MPI_Offset start[2]; + + start[0] = framei-1; + start[1] = 0; + + NCERR( ncmpi_begin_indep_data(ncid) ); + + for (int i = 0; i < n_perframe; i++) { + + if (perframe[i].type == THIS_IS_A_BIGINT) { + bigint data; + (this->*perframe[i].compute)((void*) &data); + + if (filewriter) +#if defined(LAMMPS_SMALLBIG) || defined(LAMMPS_BIGBIG) + NCERR( ncmpi_put_var1_long(ncid, perframe[i].var, start, &data) ); +#else + NCERR( ncmpi_put_var1_int(ncid, perframe[i].var, start, &data) ); +#endif + } + else { + double data; + int j = perframe[i].index; + int idim = perframe[i].dim; + + if (perframe[i].type == THIS_IS_A_COMPUTE) { + if (idim >= 0) { + modify->compute[j]->compute_vector(); + data = modify->compute[j]->vector[idim]; + } + else + data = modify->compute[j]->compute_scalar(); + } + else if (perframe[i].type == THIS_IS_A_FIX) { + if (idim >= 0) { + data = modify->fix[j]->compute_vector(idim); + } + else + data = modify->fix[j]->compute_scalar(); + } + else if (perframe[i].type == THIS_IS_A_VARIABLE) { + j = input->variable->find(perframe[i].id); + data = input->variable->compute_equal(j); + } + + if (filewriter) + NCERR( ncmpi_put_var1_double(ncid, perframe[i].var, start, &data) ); + } + } + + // write timestep header + + write_time_and_cell(); + + NCERR( ncmpi_end_indep_data(ncid) ); + + // nme = # of dump lines this proc contributes to dump + + nme = count(); + int *block_sizes = new int[comm->nprocs]; + MPI_Allgather(&nme, 1, MPI_INT, block_sizes, 1, MPI_INT, MPI_COMM_WORLD); + blocki = 0; + for (int i = 0; i < comm->me; i++) blocki += block_sizes[i]; + delete [] block_sizes; + + // insure buf is sized for packing and communicating + // use nme to insure filewriter proc can receive info from others + // limit nme*size_one to int since used as arg in MPI calls + + if (nme > maxbuf) { + if ((bigint) nme * size_one > MAXSMALLINT) + error->all(FLERR,"Too much per-proc info for dump"); + maxbuf = nme; + memory->destroy(buf); + memory->create(buf,maxbuf*size_one,"dump:buf"); + } + + // pack my data into buf + + pack(NULL); + + // each process writes its data + + write_data(nme, buf); + + // close file. this ensures data is flushed and minimizes data corruption + + closefile(); +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::write_time_and_cell() +{ + MPI_Offset start[2]; + + start[0] = framei-1; + start[1] = 0; + + MPI_Offset count[2]; + double time, cell_origin[3], cell_lengths[3], cell_angles[3]; + + time = update->ntimestep; + if (domain->triclinic == 0) { + cell_origin[0] = domain->boxlo[0]; + cell_origin[1] = domain->boxlo[1]; + cell_origin[2] = domain->boxlo[2]; + + cell_lengths[0] = domain->xprd; + cell_lengths[1] = domain->yprd; + cell_lengths[2] = domain->zprd; + + cell_angles[0] = 90; + cell_angles[1] = 90; + cell_angles[2] = 90; + } + else { + double cosalpha, cosbeta, cosgamma; + double *h = domain->h; + + cell_origin[0] = domain->boxlo[0]; + cell_origin[1] = domain->boxlo[1]; + cell_origin[2] = domain->boxlo[2]; + + cell_lengths[0] = domain->xprd; + cell_lengths[1] = sqrt(h[1]*h[1]+h[5]*h[5]); + cell_lengths[2] = sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]); + + cosalpha = (h[5]*h[4]+h[1]*h[3])/ + sqrt((h[1]*h[1]+h[5]*h[5])*(h[2]*h[2]+h[3]*h[3]+h[4]*h[4])); + cosbeta = h[4]/sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]); + cosgamma = h[5]/sqrt(h[1]*h[1]+h[5]*h[5]); + + cell_angles[0] = acos(cosalpha)*180.0/MY_PI; + cell_angles[1] = acos(cosbeta)*180.0/MY_PI; + cell_angles[2] = acos(cosgamma)*180.0/MY_PI; + } + + // Recent AMBER conventions say that nonperiodic boundaries should have + // 'cell_lengths' set to zero. + for (int dim = 0; dim < 3; dim++) { + if (!domain->periodicity[dim]) + cell_lengths[dim] = 0.0; + } + + count[0] = 1; + count[1] = 3; + if (filewriter) { + NCERR( ncmpi_put_var1_double(ncid, time_var, start, &time) ); + NCERR( ncmpi_put_vara_double(ncid, cell_origin_var, start, count, + cell_origin) ); + NCERR( ncmpi_put_vara_double(ncid, cell_lengths_var, start, count, + cell_lengths) ); + NCERR( ncmpi_put_vara_double(ncid, cell_angles_var, start, count, + cell_angles) ); + } +} + + +/* ---------------------------------------------------------------------- + write data lines to file in a block-by-block style + write head of block (mass & element name) only if has atoms of the type +------------------------------------------------------------------------- */ + +void DumpNCMPIIO::write_data(int n, double *mybuf) +{ + MPI_Offset start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + MPI_Offset stride[NC_MAX_VAR_DIMS]; + + if (!int_buffer) { + n_buffer = std::max(1, n); + int_buffer = (int *) + memory->smalloc(n_buffer*sizeof(int), "DumpNCMPIIO::int_buffer"); + double_buffer = (double *) + memory->smalloc(n_buffer*sizeof(double), "DumpNCMPIIO::double_buffer"); + } + + if (n > n_buffer) { + n_buffer = std::max(1, n); + int_buffer = (int *) + memory->srealloc(int_buffer, n_buffer*sizeof(int), + "DumpNCMPIIO::int_buffer"); + double_buffer = (double *) + memory->srealloc(double_buffer, n_buffer*sizeof(double), + "DumpNCMPIIO::double_buffer"); + } + + start[0] = framei-1; + start[1] = blocki; + start[2] = 0; + + if (n == 0) { + /* If there is no data, we need to make sure the start values don't exceed + dimension bounds. Just set them to zero. */ + start[1] = 0; + } + + count[0] = 1; + count[1] = n; + count[2] = 1; + + stride[0] = 1; + stride[1] = 1; + stride[2] = 3; + + for (int i = 0; i < n_perat; i++) { + int iaux = perat[i].field[0]; + if (iaux < 0 || iaux >= size_one) { + char errmsg[1024]; + sprintf(errmsg, "Internal error: name = %s, iaux = %i, " + "size_one = %i", perat[i].name, iaux, size_one); + error->one(FLERR,errmsg); + } + + if (vtype[iaux] == INT) { + // integers + if (perat[i].dims > 1) { + + for (int idim = 0; idim < perat[i].dims; idim++) { + iaux = perat[i].field[idim]; + + if (iaux >= 0) { + if (iaux >= size_one) { + char errmsg[1024]; + sprintf(errmsg, "Internal error: name = %s, iaux = %i, " + "size_one = %i", perat[i].name, iaux, size_one); + error->one(FLERR,errmsg); + } + + for (int j = 0; j < n; j++, iaux+=size_one) { + int_buffer[j] = mybuf[iaux]; + } + + start[2] = idim; + NCERRX( ncmpi_put_vars_int_all(ncid, perat[i].var, start, count, + stride, int_buffer), perat[i].name ); + } + } + } + else { + for (int j = 0; j < n; j++, iaux+=size_one) { + int_buffer[j] = mybuf[iaux]; + } + + NCERRX( ncmpi_put_vara_int_all(ncid, perat[i].var, start, count, + int_buffer), perat[i].name ); + } + } + else { + // doubles + if (perat[i].dims > 1) { + + for (int idim = 0; idim < perat[i].dims; idim++) { + iaux = perat[i].field[idim]; + + if (iaux >= 0) { + if (iaux >= size_one) { + char errmsg[1024]; + sprintf(errmsg, "Internal error: name = %s, iaux = %i, " + "size_one = %i", perat[i].name, iaux, size_one); + error->one(FLERR,errmsg); + } + + for (int j = 0; j < n; j++, iaux+=size_one) { + double_buffer[j] = mybuf[iaux]; + } + + start[2] = idim; + NCERRX( ncmpi_put_vars_double_all(ncid, perat[i].var, start, count, + stride, double_buffer), perat[i].name ); + } + } + } + else { + for (int j = 0; j < n; j++, iaux+=size_one) { + double_buffer[j] = mybuf[iaux]; + } + + NCERRX( ncmpi_put_vara_double_all(ncid, perat[i].var, start, count, + double_buffer), perat[i].name ); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +int DumpNCMPIIO::modify_param(int narg, char **arg) +{ + int iarg = 0; + if (strcmp(arg[iarg],"double") == 0) { + iarg++; + if (iarg >= narg) + error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword."); + if (strcmp(arg[iarg],"yes") == 0) { + double_precision = true; + } + else if (strcmp(arg[iarg],"no") == 0) { + double_precision = false; + } + else error->all(FLERR,"expected 'yes' or 'no' after 'double' keyword."); + iarg++; + return 2; + } + else if (strcmp(arg[iarg],"at") == 0) { + iarg++; + framei = force->inumeric(FLERR,arg[iarg]); + if (framei < 0) framei--; + iarg++; + return 2; + } + else if (strcmp(arg[iarg],"global") == 0) { + // "perframe" quantities, i.e. not per-atom stuff + + iarg++; + + n_perframe = narg-iarg; + perframe = new nc_perframe_t[n_perframe]; + + for (int i = 0; iarg < narg; iarg++, i++) { + int n; + char *suffix; + + if (!strcmp(arg[iarg],"step")) { + perframe[i].type = THIS_IS_A_BIGINT; + perframe[i].compute = &DumpNCMPIIO::compute_step; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strcmp(arg[iarg],"elapsed")) { + perframe[i].type = THIS_IS_A_BIGINT; + perframe[i].compute = &DumpNCMPIIO::compute_elapsed; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strcmp(arg[iarg],"elaplong")) { + perframe[i].type = THIS_IS_A_BIGINT; + perframe[i].compute = &DumpNCMPIIO::compute_elapsed_long; + strcpy(perframe[i].name, arg[iarg]); + } + else { + + n = strlen(arg[iarg]); + + if (n > 2) { + suffix = new char[n-1]; + strcpy(suffix, arg[iarg]+2); + } + else { + char errstr[1024]; + sprintf(errstr, "perframe quantity '%s' must thermo quantity or " + "compute, fix or variable", arg[iarg]); + error->all(FLERR,errstr); + } + + if (!strncmp(arg[iarg], "c_", 2)) { + int idim = -1; + char *ptr = strchr(suffix, '['); + + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all(FLERR,"Missing ']' in dump modify command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + } + + n = modify->find_compute(suffix); + if (n < 0) + error->all(FLERR,"Could not find dump modify compute ID"); + if (modify->compute[n]->peratom_flag != 0) + error->all(FLERR,"Dump modify compute ID computes per-atom info"); + if (idim >= 0 && modify->compute[n]->vector_flag == 0) + error->all(FLERR,"Dump modify compute ID does not compute vector"); + if (idim < 0 && modify->compute[n]->scalar_flag == 0) + error->all(FLERR,"Dump modify compute ID does not compute scalar"); + + perframe[i].type = THIS_IS_A_COMPUTE; + perframe[i].dim = idim; + perframe[i].index = n; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strncmp(arg[iarg], "f_", 2)) { + int idim = -1; + char *ptr = strchr(suffix, '['); + + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all(FLERR,"Missing ']' in dump modify command"); + *ptr = '\0'; + idim = ptr[1] - '1'; + } + + n = modify->find_fix(suffix); + if (n < 0) + error->all(FLERR,"Could not find dump modify fix ID"); + if (modify->fix[n]->peratom_flag != 0) + error->all(FLERR,"Dump modify fix ID computes per-atom info"); + if (idim >= 0 && modify->fix[n]->vector_flag == 0) + error->all(FLERR,"Dump modify fix ID does not compute vector"); + if (idim < 0 && modify->fix[n]->scalar_flag == 0) + error->all(FLERR,"Dump modify fix ID does not compute vector"); + + perframe[i].type = THIS_IS_A_FIX; + perframe[i].dim = idim; + perframe[i].index = n; + strcpy(perframe[i].name, arg[iarg]); + } + else if (!strncmp(arg[iarg], "v_", 2)) { + n = input->variable->find(suffix); + if (n < 0) + error->all(FLERR,"Could not find dump modify variable ID"); + if (!input->variable->equalstyle(n)) + error->all(FLERR,"Dump modify variable must be of style equal"); + + perframe[i].type = THIS_IS_A_VARIABLE; + perframe[i].dim = 1; + perframe[i].index = n; + strcpy(perframe[i].name, arg[iarg]); + strcpy(perframe[i].id, suffix); + } + else { + char errstr[1024]; + sprintf(errstr, "perframe quantity '%s' must be compute, fix or " + "variable", arg[iarg]); + error->all(FLERR,errstr); + } + + delete [] suffix; + + } + } + + return narg; + } else return 0; +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::ncerr(int err, const char *descr, int line) +{ + if (err != NC_NOERR) { + char errstr[1024]; + if (descr) { + sprintf(errstr, "NetCDF failed with error '%s' (while accessing '%s') " + " in line %i of %s.", ncmpi_strerror(err), descr, line, __FILE__); + } + else { + sprintf(errstr, "NetCDF failed with error '%s' in line %i of %s.", + ncmpi_strerror(err), line, __FILE__); + } + error->one(FLERR,errstr); + } +} + +/* ---------------------------------------------------------------------- + one method for every keyword thermo can output + called by compute() or evaluate_keyword() + compute will have already been called + set ivalue/dvalue/bivalue if value is int/double/bigint + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +void DumpNCMPIIO::compute_step(void *r) +{ + *((bigint *) r) = update->ntimestep; +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::compute_elapsed(void *r) +{ + *((bigint *) r) = update->ntimestep - update->firststep; +} + +/* ---------------------------------------------------------------------- */ + +void DumpNCMPIIO::compute_elapsed_long(void *r) +{ + *((bigint *) r) = update->ntimestep - update->beginstep; +} + +#endif /* defined(LMP_HAS_PNETCDF) */ diff --git a/src/USER-NC-DUMP/dump_nc_mpiio.h b/src/USER-NC-DUMP/dump_nc_mpiio.h new file mode 100644 index 0000000000..e7f85b4bfc --- /dev/null +++ b/src/USER-NC-DUMP/dump_nc_mpiio.h @@ -0,0 +1,140 @@ +/* ====================================================================== + LAMMPS NetCDF dump style + https://github.com/pastewka/lammps-netcdf + Lars Pastewka, lars.pastewka@kit.edu + + Copyright (2011-2013) Fraunhofer IWM + Copyright (2014) Karlsruhe Institute of Technology + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + ====================================================================== */ + +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +#if defined(LMP_HAS_PNETCDF) + +#ifdef DUMP_CLASS + +DumpStyle(nc/mpiio,DumpNCMPIIO) + +#else + +#ifndef LMP_DUMP_NC_MPIIO_H +#define LMP_DUMP_NC_MPIIO_H + +#include "dump_custom.h" + +namespace LAMMPS_NS { + +const int NC_MPIIO_FIELD_NAME_MAX = 100; +const int DUMP_NC_MPIIO_MAX_DIMS = 100; + +class DumpNCMPIIO : public DumpCustom { + public: + DumpNCMPIIO(class LAMMPS *, int, char **); + virtual ~DumpNCMPIIO(); + virtual void write(); + + private: + // per-atoms quantities (positions, velocities, etc.) + struct nc_perat_t { + int dims; // number of dimensions + int field[DUMP_NC_MPIIO_MAX_DIMS]; // field indices corresponding to the dim. + char name[NC_MPIIO_FIELD_NAME_MAX]; // field name + int var; // NetCDF variable + }; + + typedef void (DumpNCMPIIO::*funcptr_t)(void *); + + // per-frame quantities (variables, fixes or computes) + struct nc_perframe_t { + char name[NC_MPIIO_FIELD_NAME_MAX]; // field name + int var; // NetCDF variable + int type; // variable, fix, compute or callback + int index; // index in fix/compute list + funcptr_t compute; // compute function + int dim; // dimension + char id[NC_MPIIO_FIELD_NAME_MAX]; // variable id + + bigint bigint_data; // actual data + double double_data; // actual data + }; + + int framei; // current frame index + int blocki; // current block index + int ndata; // number of data blocks to expect + + bigint ntotalgr; // # of atoms + + int n_perat; // # of netcdf per-atom properties + nc_perat_t *perat; // per-atom properties + + int n_perframe; // # of global netcdf (not per-atom) fix props + nc_perframe_t *perframe; // global properties + + bool double_precision; // write everything as double precision + + bigint n_buffer; // size of buffer + int *int_buffer; // buffer for passing data to netcdf + double *double_buffer; // buffer for passing data to netcdf + + int ncid; + + int frame_dim; + int spatial_dim; + int Voigt_dim; + int atom_dim; + int cell_spatial_dim; + int cell_angular_dim; + int label_dim; + + int spatial_var; + int cell_spatial_var; + int cell_angular_var; + + int time_var; + int cell_origin_var; + int cell_lengths_var; + int cell_angles_var; + + virtual void openfile(); + void closefile(); + void write_time_and_cell(); + virtual void write_data(int, double *); + void write_prmtop(); + + virtual int modify_param(int, char **); + + void ncerr(int, const char *, int); + + void compute_step(void *); + void compute_elapsed(void *); + void compute_elapsed_long(void *); +}; + +} + +#endif +#endif +#endif /* defined(LMP_HAS_PNETCDF) */