ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5458 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2011-01-04 19:13:28 +00:00
parent d5aa03447b
commit fceac703e5
1 changed files with 366 additions and 366 deletions
Download Patch File Download Diff File Expand all files Collapse all files

732
src/dump_dcd.cpp
View File

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