Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/USER-NETCDF/dump_netcdf.cpp
4beccf508f9clammm-devel

Last Change
Raw File

dump_netcdf.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Lars Pastewka (University of Freiburg)
------------------------------------------------------------------------- */
#if defined(LMP_HAS_NETCDF)
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <netcdf.h>
#include "dump_netcdf.h"
#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 "update.h"
#include "universe.h"
#include "variable.h"
#include "force.h"
#include "output.h"
#include "thermo.h"
using namespace LAMMPS_NS;
using namespace MathConst;
enum{INT,FLOAT,BIGINT}; // same as in thermo.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__)
#if !defined(NC_64BIT_DATA)
#define NC_64BIT_DATA NC_64BIT_OFFSET
#endif
/* ---------------------------------------------------------------------- */
DumpNetCDF::DumpNetCDF(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_buffer = 0;
int_buffer = NULL;
double_buffer = NULL;
double_precision = false;
thermo = false;
thermovar = NULL;
framei = 0;
}
/* ---------------------------------------------------------------------- */
DumpNetCDF::~DumpNetCDF()
{
closefile();
delete [] perat;
if (thermovar) delete [] thermovar;
if (int_buffer) memory->sfree(int_buffer);
if (double_buffer) memory->sfree(double_buffer);
}
/* ---------------------------------------------------------------------- */
void DumpNetCDF::openfile()
{
if (thermo && !singlefile_opened) {
if (thermovar) delete [] thermovar;
thermovar = new int[output->thermo->nfield];
}
// 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.
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
if (thermo) {
Thermo *th = output->thermo;
for (int i = 0; i < th->nfield; i++) {
NCERRX( nc_inq_varid(ncid, th->keyword[i], &thermovar[i]),
th->keyword[i] );
}
}
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_DATA, &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
if (thermo) {
Thermo *th = output->thermo;
for (int i = 0; i < th->nfield; i++) {
if (th->vtype[i] == FLOAT) {
NCERRX( nc_def_var(ncid, th->keyword[i], NC_DOUBLE, 1, dims,
&thermovar[i]), th->keyword[i] );
}
else if (th->vtype[i] == INT) {
NCERRX( nc_def_var(ncid, th->keyword[i], NC_INT, 1, dims,
&thermovar[i]), th->keyword[i] );
}
else if (th->vtype[i] == BIGINT) {
NCERRX( nc_def_var(ncid, th->keyword[i], NC_LONG, 1, dims,
&thermovar[i]), th->keyword[i] );
}
}
}
// 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 DumpNetCDF::closefile()
{
if (filewriter && singlefile_opened) {
NCERR( nc_close(ncid) );
singlefile_opened = 0;
// append next time DumpNetCDF::openfile is called
append_flag = 1;
// write to next frame upon next open
framei++;
}
}
/* ---------------------------------------------------------------------- */
void DumpNetCDF::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;
if (thermo) {
Thermo *th = output->thermo;
for (int i = 0; i < th->nfield; i++) {
th->call_vfunc(i);
if (filewriter) {
if (th->vtype[i] == FLOAT) {
NCERRX( nc_put_var1_double(ncid, thermovar[i], start,
&th->dvalue),
th->keyword[i] );
}
else if (th->vtype[i] == INT) {
NCERRX( nc_put_var1_int(ncid, thermovar[i], start, &th->ivalue),
th->keyword[i] );
}
else if (th->vtype[i] == BIGINT) {
#if defined(LAMMPS_SMALLBIG) || defined(LAMMPS_BIGBIG)
NCERRX( nc_put_var1_long(ncid, thermovar[i], start, &th->bivalue),
th->keyword[i] );
#else
NCERRX( nc_put_var1_int(ncid, thermovar[i], start, &th->bivalue),
th->keyword[i] );
#endif
}
}
}
}
// call write of superclass
Dump::write();
// close file. this ensures data is flushed and mimized data corruption
closefile();
}
/* ---------------------------------------------------------------------- */
void DumpNetCDF::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 DumpNetCDF::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),"dump::int_buffer");
double_buffer = (double *)
memory->smalloc(n*sizeof(double),"dump::double_buffer");
}
if (n > n_buffer) {
n_buffer = n;
int_buffer = (int *)
memory->srealloc(int_buffer, n*sizeof(int),"dump::int_buffer");
double_buffer = (double *)
memory->srealloc(double_buffer, n*sizeof(double),"dump::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 DumpNetCDF::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],"thermo") == 0) {
iarg++;
if (iarg >= narg)
error->all(FLERR,"expected 'yes' or 'no' after 'thermo' keyword.");
if (strcmp(arg[iarg],"yes") == 0) {
thermo = true;
}
else if (strcmp(arg[iarg],"no") == 0) {
thermo = false;
}
else error->all(FLERR,"expected 'yes' or 'no' after 'thermo' keyword.");
iarg++;
return 2;
} else return 0;
}
/* ---------------------------------------------------------------------- */
void DumpNetCDF::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 DumpNetCDF::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);
}
}
#endif /* defined(LMP_HAS_NETCDF) */
c4science · Help