diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp index 0ae2b707ea..6794cfbc64 100644 --- a/src/dump_dcd.cpp +++ b/src/dump_dcd.cpp @@ -1,366 +1,366 @@ -/* ---------------------------------------------------------------------- - 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: Naveen Michaud-Agrawal (Johns Hopkins U) - Axel Kohlmeyer (Temple U), support for groups -------------------------------------------------------------------------- */ - -#include "math.h" -#include "inttypes.h" -#include "stdio.h" -#include "time.h" -#include "string.h" -#include "dump_dcd.h" -#include "domain.h" -#include "atom.h" -#include "update.h" -#include "output.h" -#include "group.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define NFILE_POS 8L -#define NSTEP_POS 20L - -// necessary to set SEEK params b/c MPI-2 messes with these settings - -#ifndef SEEK_SET -#define SEEK_SET 0 -#define SEEK_CUR 1 -#define SEEK_END 2 -#endif - -/* ---------------------------------------------------------------------- */ - -static inline void fwrite_int32(FILE* fd, uint32_t i) -{ - fwrite(&i,sizeof(uint32_t),1,fd); -} - -/* ---------------------------------------------------------------------- */ - -DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) -{ - if (narg != 5) error->all("Illegal dump dcd command"); - if (binary || compressed || multifile || multiproc) - error->all("Invalid dump dcd filename"); - - size_one = 3; - sort_flag = 1; - sortcol = 0; - - unwrap_flag = 0; - format_default = NULL; - - // allocate global array for atom coords - - if (igroup == 0) natoms = static_cast (atom->natoms); - else natoms = static_cast (group->count(igroup)); - if (natoms <= 0) error->all("Invalid natoms for dump dcd"); - - coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); - xf = &coords[0*natoms]; - yf = &coords[1*natoms]; - zf = &coords[2*natoms]; - - openfile(); - headerflag = 0; - nevery_save = 0; - ntotal = 0; -} - -/* ---------------------------------------------------------------------- */ - -DumpDCD::~DumpDCD() -{ - memory->sfree(coords); -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::init_style() -{ - if (sort_flag == 0 || sortcol != 0) - error->all("Dump dcd requires sorting by atom ID"); - - // check that dump frequency has not changed and is not a variable - - int idump; - for (idump = 0; idump < output->ndump; idump++) - if (strcmp(id,output->dump[idump]->id) == 0) break; - if (output->every_dump[idump] == 0) - error->all("Cannot use variable every setting for dump dcd"); - - if (nevery_save == 0) nevery_save = output->every_dump[idump]; - else if (nevery_save != output->every_dump[idump]) - error->all("Cannot change dump_modify every for dump dcd"); -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::openfile() -{ - if (me == 0) { - fp = fopen(filename,"wb"); - if (fp == NULL) error->one("Cannot open dump file"); - } -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::write_header(int n) -{ - if (n != natoms) error->all("Dump dcd of non-matching # of atoms"); - - // first time, write header for entire file - - if (headerflag == 0) { - if (me == 0) write_dcd_header("Written by LAMMPS"); - headerflag = 1; - nframes = 0; - } - - // dim[] = size and angle cosines of orthogonal or triclinic box - // dim[0] = a = length of unit cell vector along x-axis - // dim[1] = gamma = cosine of angle between a and b - // dim[2] = b = length of unit cell vector in xy-plane - // dim[3] = beta = cosine of angle between a and c - // dim[4] = alpha = cosine of angle between b and c - // dim[5] = c = length of final unit cell vector - // 48 = 6 doubles - - double dim[6]; - if (domain->triclinic) { - double *h = domain->h; - double alen = h[0]; - double blen = sqrt(h[5]*h[5] + h[1]*h[1]); - double clen = sqrt(h[4]*h[4] + h[3]*h[3] + h[2]*h[2]); - dim[0] = alen; - dim[2] = blen; - dim[5] = clen; - dim[4] = (h[5]*h[4] + h[1]*h[3]) / blen/clen; // alpha - dim[3] = (h[0]*h[4]) / alen/clen; // beta - dim[1] = (h[0]*h[5]) / alen/blen; // gamma - } else { - dim[0] = domain->xprd; - dim[2] = domain->yprd; - dim[5] = domain->zprd; - dim[1] = dim[3] = dim[4] = 0.0; - } - - if (me == 0) { - uint32_t out_integer = 48; - fwrite_int32(fp,out_integer); - fwrite(dim,out_integer,1,fp); - fwrite_int32(fp,out_integer); - if (flush_flag) fflush(fp); - } -} - -/* ---------------------------------------------------------------------- */ - -int DumpDCD::count() -{ - if (igroup == 0) return atom->nlocal; - - int *mask = atom->mask; - int nlocal = atom->nlocal; - - int m = 0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) m++; - return m; -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::pack(int *ids) -{ - int m,n; - - int *tag = atom->tag; - double **x = atom->x; - int *image = atom->image; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - m = n = 0; - if (unwrap_flag) { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double xy = domain->xy; - double xz = domain->xz; - double yz = domain->yz; - - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - int ix = (image[i] & 1023) - 512; - int iy = (image[i] >> 10 & 1023) - 512; - int iz = (image[i] >> 20) - 512; - - if (domain->triclinic) { - buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz; - buf[m++] = x[i][1] + iy * yprd + iz * yz; - buf[m++] = x[i][2] + iz * zprd; - } else { - buf[m++] = x[i][0] + ix * xprd; - buf[m++] = x[i][1] + iy * yprd; - buf[m++] = x[i][2] + iz * zprd; - } - ids[n++] = tag[i]; - } - } - - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - ids[n++] = tag[i]; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::write_data(int n, double *mybuf) -{ - // copy buf atom coords into 3 global arrays - - int m = 0; - for (int i = 0; i < n; i++) { - xf[ntotal] = mybuf[m++]; - yf[ntotal] = mybuf[m++]; - zf[ntotal] = mybuf[m++]; - ntotal++; - } - - // if last chunk of atoms in this snapshot, write global arrays to file - - if (ntotal == natoms) { - write_frame(); - ntotal = 0; - } -} - -/* ---------------------------------------------------------------------- */ - -int DumpDCD::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"unwrap") == 0) { - if (narg < 2) error->all("Illegal dump_modify command"); - if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; - else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; - else error->all("Illegal dump_modify command"); - return 2; - } - return 0; -} - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory in buf and global coords array -------------------------------------------------------------------------- */ - -double DumpDCD::memory_usage() -{ - double bytes = Dump::memory_usage(); - bytes += 3*natoms * sizeof(float); - return bytes; -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::write_frame() -{ - // write coords - - uint32_t out_integer = natoms*sizeof(float); - fwrite_int32(fp,out_integer); - fwrite(xf,out_integer,1,fp); - fwrite_int32(fp,out_integer); - fwrite_int32(fp,out_integer); - fwrite(yf,out_integer,1,fp); - fwrite_int32(fp,out_integer); - fwrite_int32(fp,out_integer); - fwrite(zf,out_integer,1,fp); - fwrite_int32(fp,out_integer); - - // update NFILE and NSTEP fields in DCD header - - nframes++; - out_integer = nframes; - fseek(fp,NFILE_POS,SEEK_SET); - fwrite_int32(fp,out_integer); - out_integer = update->ntimestep; - fseek(fp,NSTEP_POS,SEEK_SET); - fwrite_int32(fp,out_integer); - fseek(fp,0,SEEK_END); -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::write_dcd_header(const char *remarks) -{ - uint32_t out_integer; - float out_float; - char title_string[200]; - time_t cur_time; - struct tm *tmbuf; - - out_integer = 84; - fwrite_int32(fp,out_integer); - strcpy(title_string,"CORD"); - fwrite(title_string,4,1,fp); - fwrite_int32(fp,0); // NFILE = # of snapshots in file - fwrite_int32(fp,update->ntimestep); // START = timestep of first snapshot - fwrite_int32(fp,nevery_save); // SKIP = interval between snapshots - fwrite_int32(fp,update->ntimestep); // NSTEP = timestep of last snapshot - fwrite_int32(fp,0); // NAMD writes NSTEP or ISTART - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - out_float = update->dt; - fwrite(&out_float,sizeof(float),1,fp); - fwrite_int32(fp,1); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,0); - fwrite_int32(fp,24); // pretend to be Charmm version 24 - fwrite_int32(fp,84); - fwrite_int32(fp,164); - fwrite_int32(fp,2); - strncpy(title_string,remarks,80); - title_string[79] = '\0'; - fwrite(title_string,80,1,fp); - cur_time=time(NULL); - tmbuf=localtime(&cur_time); - memset(title_string,' ',81); - strftime(title_string,80,"REMARKS Created %d %B,%Y at %R",tmbuf); - fwrite(title_string,80,1,fp); - fwrite_int32(fp,164); - fwrite_int32(fp,4); - fwrite_int32(fp,natoms); // number of atoms in each snapshot - fwrite_int32(fp,4); - if (flush_flag) fflush(fp); -} +/* ---------------------------------------------------------------------- + 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: Naveen Michaud-Agrawal (Johns Hopkins U) + Axel Kohlmeyer (Temple U), support for groups +------------------------------------------------------------------------- */ + +#include "math.h" +#include "inttypes.h" +#include "stdio.h" +#include "time.h" +#include "string.h" +#include "dump_dcd.h" +#include "domain.h" +#include "atom.h" +#include "update.h" +#include "output.h" +#include "group.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define NFILE_POS 8L +#define NSTEP_POS 20L + +// necessary to set SEEK params b/c MPI-2 messes with these settings + +#ifndef SEEK_SET +#define SEEK_SET 0 +#define SEEK_CUR 1 +#define SEEK_END 2 +#endif + +/* ---------------------------------------------------------------------- */ + +static inline void fwrite_int32(FILE* fd, uint32_t i) +{ + fwrite(&i,sizeof(uint32_t),1,fd); +} + +/* ---------------------------------------------------------------------- */ + +DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) +{ + if (narg != 5) error->all("Illegal dump dcd command"); + if (binary || compressed || multifile || multiproc) + error->all("Invalid dump dcd filename"); + + size_one = 3; + sort_flag = 1; + sortcol = 0; + + unwrap_flag = 0; + format_default = NULL; + + // allocate global array for atom coords + + if (igroup == 0) natoms = static_cast (atom->natoms); + else natoms = static_cast (group->count(igroup)); + if (natoms <= 0) error->all("Invalid natoms for dump dcd"); + + coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); + xf = &coords[0*natoms]; + yf = &coords[1*natoms]; + zf = &coords[2*natoms]; + + openfile(); + headerflag = 0; + nevery_save = 0; + ntotal = 0; +} + +/* ---------------------------------------------------------------------- */ + +DumpDCD::~DumpDCD() +{ + memory->sfree(coords); +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::init_style() +{ + if (sort_flag == 0 || sortcol != 0) + error->all("Dump dcd requires sorting by atom ID"); + + // check that dump frequency has not changed and is not a variable + + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + if (output->every_dump[idump] == 0) + error->all("Cannot use variable every setting for dump dcd"); + + if (nevery_save == 0) nevery_save = output->every_dump[idump]; + else if (nevery_save != output->every_dump[idump]) + error->all("Cannot change dump_modify every for dump dcd"); +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::openfile() +{ + if (me == 0) { + fp = fopen(filename,"wb"); + if (fp == NULL) error->one("Cannot open dump file"); + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::write_header(int n) +{ + if (n != natoms) error->all("Dump dcd of non-matching # of atoms"); + + // first time, write header for entire file + + if (headerflag == 0) { + if (me == 0) write_dcd_header("Written by LAMMPS"); + headerflag = 1; + nframes = 0; + } + + // dim[] = size and angle cosines of orthogonal or triclinic box + // dim[0] = a = length of unit cell vector along x-axis + // dim[1] = gamma = cosine of angle between a and b + // dim[2] = b = length of unit cell vector in xy-plane + // dim[3] = beta = cosine of angle between a and c + // dim[4] = alpha = cosine of angle between b and c + // dim[5] = c = length of final unit cell vector + // 48 = 6 doubles + + double dim[6]; + if (domain->triclinic) { + double *h = domain->h; + double alen = h[0]; + double blen = sqrt(h[5]*h[5] + h[1]*h[1]); + double clen = sqrt(h[4]*h[4] + h[3]*h[3] + h[2]*h[2]); + dim[0] = alen; + dim[2] = blen; + dim[5] = clen; + dim[4] = (h[5]*h[4] + h[1]*h[3]) / blen/clen; // alpha + dim[3] = (h[0]*h[4]) / alen/clen; // beta + dim[1] = (h[0]*h[5]) / alen/blen; // gamma + } else { + dim[0] = domain->xprd; + dim[2] = domain->yprd; + dim[5] = domain->zprd; + dim[1] = dim[3] = dim[4] = 0.0; + } + + if (me == 0) { + uint32_t out_integer = 48; + fwrite_int32(fp,out_integer); + fwrite(dim,out_integer,1,fp); + fwrite_int32(fp,out_integer); + if (flush_flag) fflush(fp); + } +} + +/* ---------------------------------------------------------------------- */ + +int DumpDCD::count() +{ + if (igroup == 0) return atom->nlocal; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int m = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) m++; + return m; +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::pack(int *ids) +{ + int m,n; + + int *tag = atom->tag; + double **x = atom->x; + int *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + m = n = 0; + if (unwrap_flag) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xy = domain->xy; + double xz = domain->xz; + double yz = domain->yz; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + int ix = (image[i] & 1023) - 512; + int iy = (image[i] >> 10 & 1023) - 512; + int iz = (image[i] >> 20) - 512; + + if (domain->triclinic) { + buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz; + buf[m++] = x[i][1] + iy * yprd + iz * yz; + buf[m++] = x[i][2] + iz * zprd; + } else { + buf[m++] = x[i][0] + ix * xprd; + buf[m++] = x[i][1] + iy * yprd; + buf[m++] = x[i][2] + iz * zprd; + } + ids[n++] = tag[i]; + } + } + + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + ids[n++] = tag[i]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::write_data(int n, double *mybuf) +{ + // copy buf atom coords into 3 global arrays + + int m = 0; + for (int i = 0; i < n; i++) { + xf[ntotal] = mybuf[m++]; + yf[ntotal] = mybuf[m++]; + zf[ntotal] = mybuf[m++]; + ntotal++; + } + + // if last chunk of atoms in this snapshot, write global arrays to file + + if (ntotal == natoms) { + write_frame(); + ntotal = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +int DumpDCD::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"unwrap") == 0) { + if (narg < 2) error->all("Illegal dump_modify command"); + if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; + else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; + else error->all("Illegal dump_modify command"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory in buf and global coords array +------------------------------------------------------------------------- */ + +double DumpDCD::memory_usage() +{ + double bytes = Dump::memory_usage(); + bytes += 3*natoms * sizeof(float); + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::write_frame() +{ + // write coords + + uint32_t out_integer = natoms*sizeof(float); + fwrite_int32(fp,out_integer); + fwrite(xf,out_integer,1,fp); + fwrite_int32(fp,out_integer); + fwrite_int32(fp,out_integer); + fwrite(yf,out_integer,1,fp); + fwrite_int32(fp,out_integer); + fwrite_int32(fp,out_integer); + fwrite(zf,out_integer,1,fp); + fwrite_int32(fp,out_integer); + + // update NFILE and NSTEP fields in DCD header + + nframes++; + out_integer = nframes; + fseek(fp,NFILE_POS,SEEK_SET); + fwrite_int32(fp,out_integer); + out_integer = update->ntimestep; + fseek(fp,NSTEP_POS,SEEK_SET); + fwrite_int32(fp,out_integer); + fseek(fp,0,SEEK_END); +} + +/* ---------------------------------------------------------------------- */ + +void DumpDCD::write_dcd_header(const char *remarks) +{ + uint32_t out_integer; + float out_float; + char title_string[200]; + time_t cur_time; + struct tm *tmbuf; + + out_integer = 84; + fwrite_int32(fp,out_integer); + strcpy(title_string,"CORD"); + fwrite(title_string,4,1,fp); + fwrite_int32(fp,0); // NFILE = # of snapshots in file + fwrite_int32(fp,update->ntimestep); // START = timestep of first snapshot + fwrite_int32(fp,nevery_save); // SKIP = interval between snapshots + fwrite_int32(fp,update->ntimestep); // NSTEP = timestep of last snapshot + fwrite_int32(fp,0); // NAMD writes NSTEP or ISTART + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + out_float = update->dt; + fwrite(&out_float,sizeof(float),1,fp); + fwrite_int32(fp,1); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,0); + fwrite_int32(fp,24); // pretend to be Charmm version 24 + fwrite_int32(fp,84); + fwrite_int32(fp,164); + fwrite_int32(fp,2); + strncpy(title_string,remarks,80); + title_string[79] = '\0'; + fwrite(title_string,80,1,fp); + cur_time=time(NULL); + tmbuf=localtime(&cur_time); + memset(title_string,' ',81); + strftime(title_string,80,"REMARKS Created %d %B,%Y at %R",tmbuf); + fwrite(title_string,80,1,fp); + fwrite_int32(fp,164); + fwrite_int32(fp,4); + fwrite_int32(fp,natoms); // number of atoms in each snapshot + fwrite_int32(fp,4); + if (flush_flag) fflush(fp); +}