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@12773 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2014-11-24 17:15:05 +00:00
parent a52ad2d6dd
commit 162e8050ad
10 changed files with 2024 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/Section_commands.txt
View File

@ -583,7 +583,9 @@ package"_Section_start.html#start_3.
"atc"_fix_atc.html,
"ave/spatial/sphere"_fix_ave_spatial_sphere.html,
"colvars"_fix_colvars.html,
"gle"_fix_gle.html,
"imd"_fix_imd.html,
"ipi"_fix_ipi.html,
"langevin/eff"_fix_langevin_eff.html,
"lb/fluid"_fix_lb_fluid.html,
"lb/momentum"_fix_lb_momentum.html,

1
doc/Section_intro.txt
View File

@ -254,6 +254,7 @@ molecular dynamics options:
"atom-to-continuum coupling"_fix_atc.html with finite elements
coupled rigid body integration via the "POEMS"_fix_poems.html library
"QM/MM coupling"_fix_qmmm.html
"path-integral molecular dynamics (PIMD)"_fix_ipi.html
"grand canonical Monte Carlo"_fix_gcmc.html insertions/deletions
"Direct Simulation Monte Carlo"_pair_dsmc.html for low-density fluids
"Peridynamics mesoscale modeling"_pair_peri.html

1
src/MPIIO/dump_atom_mpiio.cpp
View File

@ -34,7 +34,6 @@ using namespace LAMMPS_NS;
#define DUMP_BUF_CHUNK_SIZE 16384
#define DUMP_BUF_INCREMENT_SIZE 4096
/* ---------------------------------------------------------------------- */
DumpAtomMPIIO::DumpAtomMPIIO(LAMMPS *lmp, int narg, char **arg) :

512
src/MPIIO/dump_cfg_mpiio.cpp Normal file
View File

@ -0,0 +1,512 @@
/* ----------------------------------------------------------------------
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: Paul Coffman (IBM)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "dump_cfg_mpiio.h"
#include "atom.h"
#include "domain.h"
#include "comm.h"
#include "modify.h"
#include "compute.h"
#include "input.h"
#include "fix.h"
#include "variable.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#ifdef LMP_USER_IO_TIMER
#include <sys/times.h>
#include <hwi/include/bqc/A2_inlines.h>
#include <stdlib.h>
long dumpCFGTimestamps[10];
#endif
#if defined(_OPENMP)
#include <omp.h>
#endif
using namespace LAMMPS_NS;
#define MAX_TEXT_HEADER_SIZE 4096
#define DUMP_BUF_CHUNK_SIZE 16384
#define DUMP_BUF_INCREMENT_SIZE 4096
enum{INT,DOUBLE,STRING,BIGINT}; // same as in DumpCustom
#define UNWRAPEXPAND 10.0
#define ONEFIELD 32
#define DELTA 1048576
/* ---------------------------------------------------------------------- */
DumpCFGMPIIO::DumpCFGMPIIO(LAMMPS *lmp, int narg, char **arg) :
DumpCFG(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
DumpCFGMPIIO::~DumpCFGMPIIO()
{
if (multifile == 0) MPI_File_close(&mpifh);
}
/* ---------------------------------------------------------------------- */
void DumpCFGMPIIO::openfile()
{
if (singlefile_opened) { // single file already opened, so just return after resetting filesize
mpifo = currentFileSize;
MPI_File_set_size(mpifh,mpifo+headerSize+sumFileSize);
currentFileSize = mpifo+headerSize+sumFileSize;
return;
}
if (multifile == 0) singlefile_opened = 1;
// if one file per timestep, replace '*' with current timestep
filecurrent = filename;
if (multifile) {
char *filestar = filecurrent;
filecurrent = new char[strlen(filestar) + 16];
char *ptr = strchr(filestar,'*');
*ptr = '\0';
if (padflag == 0)
sprintf(filecurrent,"%s" BIGINT_FORMAT "%s",
filestar,update->ntimestep,ptr+1);
else {
char bif[8],pad[16];
strcpy(bif,BIGINT_FORMAT);
sprintf(pad,"%%s%%0%d%s%%s",padflag,&bif[1]);
sprintf(filecurrent,pad,filestar,update->ntimestep,ptr+1);
}
*ptr = '*';
}
if (append_flag) { // append open
int err = MPI_File_open( world, filecurrent, MPI_MODE_CREATE | MPI_MODE_APPEND | MPI_MODE_WRONLY , MPI_INFO_NULL, &mpifh);
if (err != MPI_SUCCESS) {
char str[128];
sprintf(str,"Cannot open dump file %s",filecurrent);
error->one(FLERR,str);
}
int myrank;
MPI_Comm_rank(world,&myrank);
if (myrank == 0)
MPI_File_get_size(mpifh,&mpifo);
MPI_Bcast(&mpifo, 1, MPI_LMP_BIGINT, 0, world);
MPI_File_set_size(mpifh,mpifo+headerSize+sumFileSize);
currentFileSize = mpifo+headerSize+sumFileSize;
}
else { // replace open
int err = MPI_File_open( world, filecurrent, MPI_MODE_CREATE | MPI_MODE_WRONLY , MPI_INFO_NULL, &mpifh);
if (err != MPI_SUCCESS) {
char str[128];
sprintf(str,"Cannot open dump file %s",filecurrent);
error->one(FLERR,str);
}
mpifo = 0;
MPI_File_set_size(mpifh,(MPI_Offset) (headerSize+sumFileSize));
currentFileSize = (headerSize+sumFileSize);
}
}
/* ---------------------------------------------------------------------- */
void DumpCFGMPIIO::write()
{
#ifdef LMP_USER_IO_TIMER
long startTimeBase, endTimeBase;
MPI_Barrier(world); // timestamp barrier
if (me == 0)
startTimeBase = GetTimeBase();
#endif
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} else {
boxxlo = domain->boxlo_bound[0];
boxxhi = domain->boxhi_bound[0];
boxylo = domain->boxlo_bound[1];
boxyhi = domain->boxhi_bound[1];
boxzlo = domain->boxlo_bound[2];
boxzhi = domain->boxhi_bound[2];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
}
// nme = # of dump lines this proc contributes to dump
nme = count();
// ntotal = total # of dump lines in snapshot
// nmax = max # of dump lines on any proc
bigint bnme = nme;
MPI_Allreduce(&bnme,&ntotal,1,MPI_LMP_BIGINT,MPI_SUM,world);
int nmax;
MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world);
// write timestep header
// for multiproc,
// nheader = # of lines in this file via Allreduce on clustercomm
bigint nheader = ntotal;
// insure filewriter proc can receive everyone's info
// limit nmax*size_one to int since used as arg in MPI_Rsend() below
// pack my data into buf
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nmax > maxbuf) {
if ((bigint) nmax * size_one > MAXSMALLINT)
error->all(FLERR,"Too much per-proc info for dump");
maxbuf = nmax;
memory->destroy(buf);
memory->create(buf,(maxbuf*size_one),"dump:buf");
}
if (sort_flag && sortcol == 0 && nmax > maxids) {
maxids = nmax;
memory->destroy(ids);
memory->create(ids,maxids,"dump:ids");
}
if (sort_flag && sortcol == 0) pack(ids);
else pack(NULL);
if (sort_flag) sort();
// determine how much data needs to be written for setting the file size and prepocess it prior to writing
performEstimate = 1;
write_header(nheader);
write_data(nme,buf);
MPI_Bcast(&sumFileSize, 1, MPI_LMP_BIGINT, (nprocs-1), world);
#ifdef LMP_USER_IO_TIMER
MPI_Barrier(world); // timestamp barrier
dumpCFGTimestamps[0] = GetTimeBase();
#endif
openfile();
#ifdef LMP_USER_IO_TIMER
MPI_Barrier(world); // timestamp barrier
dumpCFGTimestamps[1] = GetTimeBase();
#endif
performEstimate = 0;
write_header(nheader); // mpifo now points to end of header info
#ifdef LMP_USER_IO_TIMER
MPI_Barrier(world); // timestamp barrier
dumpCFGTimestamps[2] = GetTimeBase();
#endif
// now actually write the data
performEstimate = 0;
write_data(nme,buf);
#ifdef LMP_USER_IO_TIMER
MPI_Barrier(world); // timestamp barrier
dumpCFGTimestamps[3] = GetTimeBase();
#endif
if (multifile) MPI_File_close(&mpifh);
if (multifile) delete [] filecurrent;
#ifdef LMP_USER_IO_TIMER
MPI_Barrier(world); // timestamp barrier
dumpCFGTimestamps[4] = GetTimeBase();
if (me == 0) {
endTimeBase = GetTimeBase();
printf("total dump cycles: %ld - estimates and setup: %ld openfile: %ld write header: %ld write data: %ld close file: %ld\n",(long) (endTimeBase-startTimeBase),(long) (dumpCFGTimestamps[0]-startTimeBase),(long) (dumpCFGTimestamps[1]-dumpCFGTimestamps[0]),(long) (dumpCFGTimestamps[2]-dumpCFGTimestamps[1]),(long) (dumpCFGTimestamps[3]-dumpCFGTimestamps[2]),(long) (dumpCFGTimestamps[4]-dumpCFGTimestamps[3]));
}
#endif
}
/* ---------------------------------------------------------------------- */
void DumpCFGMPIIO::init_style()
{
if (multifile == 0 && !multifile_override)
error->all(FLERR,"Dump cfg requires one snapshot per file");
DumpCustom::init_style();
// setup function ptrs
write_choice = &DumpCFGMPIIO::write_string;
}
/* ---------------------------------------------------------------------- */
void DumpCFGMPIIO::write_header(bigint n)
{
// set scale factor used by AtomEye for CFG viz
// default = 1.0
// for peridynamics, set to pre-computed PD scale factor
// so PD particles mimic C atoms
// for unwrapped coords, set to UNWRAPEXPAND (10.0)
// so molecules are not split across periodic box boundaries
MPI_Status mpiStatus;
if (performEstimate) {
headerBuffer = (char *) malloc(MAX_TEXT_HEADER_SIZE);
headerSize = 0;
double scale = 1.0;
if (atom->peri_flag) scale = atom->pdscale;
else if (unwrapflag == 1) scale = UNWRAPEXPAND;
char str[64];
sprintf(str,"Number of particles = %s\n",BIGINT_FORMAT);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),str,n);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"A = %g Angstrom (basic length-scale)\n",scale);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(1,1) = %g A\n",domain->xprd);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(1,2) = 0 A \n");
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(1,3) = 0 A \n");
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(2,1) = %g A \n",domain->xy);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(2,2) = %g A\n",domain->yprd);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(2,3) = 0 A \n");
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(3,1) = %g A \n",domain->xz);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(3,2) = %g A \n",domain->yz);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"H0(3,3) = %g A\n",domain->zprd);
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),".NO_VELOCITY.\n");
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"entry_count = %d\n",nfield-2);
for (int i = 0; i < nfield-5; i++)
headerSize += sprintf(((char*)&((char*)headerBuffer)[headerSize]),"auxiliary[%d] = %s\n",i,auxname[i]);
}
else { // write data
if (me == 0)
MPI_File_write_at(mpifh,mpifo,headerBuffer,headerSize,MPI_CHAR,&mpiStatus);
mpifo += headerSize;
free(headerBuffer);
}
}
#if defined(_OPENMP)
/* ----------------------------------------------------------------------
convert mybuf of doubles to one big formatted string in sbuf
return -1 if strlen exceeds an int, since used as arg in MPI calls in Dump
------------------------------------------------------------------------- */
int DumpCFGMPIIO::convert_string_omp(int n, double *mybuf)
{
MPI_Status mpiStatus;
char **mpifh_buffer_line_per_thread;
int mpifhStringCount;
int *mpifhStringCountPerThread, *bufOffset, *bufRange, *bufLength;
mpifhStringCount = 0;
int nthreads = omp_get_max_threads();
if (nthreads > n) { // call serial version
convert_string(n,mybuf);
}
else {
memory->create(mpifhStringCountPerThread,nthreads,"dump:mpifhStringCountPerThread");
mpifh_buffer_line_per_thread = (char **) malloc(nthreads*sizeof(char*));
memory->create(bufOffset,nthreads,"dump:bufOffset");
memory->create(bufRange,nthreads,"dump:bufRange");
memory->create(bufLength,nthreads,"dump:bufLength");
int i=0;
for (i=0;i<(nthreads-1);i++) {
mpifhStringCountPerThread[i] = 0;
bufOffset[i] = (int) (i*(int)(floor((double)n/(double)nthreads))*size_one);
bufRange[i] = (int)(floor((double)n/(double)nthreads));
bufLength[i] = DUMP_BUF_CHUNK_SIZE;
mpifh_buffer_line_per_thread[i] = (char *) malloc(DUMP_BUF_CHUNK_SIZE * sizeof(char));
mpifh_buffer_line_per_thread[i][0] = '\0';
}
mpifhStringCountPerThread[i] = 0;
bufOffset[i] = (int) (i*(int)(floor((double)n/(double)nthreads))*size_one);
bufRange[i] = n-(i*(int)(floor((double)n/(double)nthreads)));
bufLength[i] = DUMP_BUF_CHUNK_SIZE;
mpifh_buffer_line_per_thread[i] = (char *) malloc(DUMP_BUF_CHUNK_SIZE * sizeof(char));
mpifh_buffer_line_per_thread[i][0] = '\0';
#pragma omp parallel default(none) shared(bufOffset, bufRange, bufLength, mpifhStringCountPerThread, mpifh_buffer_line_per_thread, mybuf)
{
int tid = omp_get_thread_num();
int m=0;
if (unwrapflag == 0) {
for (int i = 0; i < bufRange[tid]; i++) {
if ((bufLength[tid] - mpifhStringCountPerThread[tid]) < DUMP_BUF_INCREMENT_SIZE) {
mpifh_buffer_line_per_thread[tid] = (char *) realloc(mpifh_buffer_line_per_thread[tid],(mpifhStringCountPerThread[tid]+DUMP_BUF_CHUNK_SIZE) * sizeof(char));
bufLength[tid] = (mpifhStringCountPerThread[tid]+DUMP_BUF_CHUNK_SIZE) * sizeof(char);
}
for (int j = 0; j < size_one; j++) {
if (j == 0) {
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),"%f \n",(mybuf[bufOffset[tid]+m]));
} else if (j == 1) {
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),"%s \n",typenames[(int) mybuf[bufOffset[tid]+m]]);
} else if (j >= 2) {
if (vtype[j] == INT)
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],static_cast<int> (mybuf[bufOffset[tid]+m]));
else if (vtype[j] == DOUBLE)
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],mybuf[bufOffset[tid]+m]);
else if (vtype[j] == STRING)
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],typenames[(int) mybuf[bufOffset[tid]+m]]);
else if (vtype[j] == BIGINT)
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],static_cast<bigint> (mybuf[bufOffset[tid]+m]));
}
m++;
} // for j
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),"\n");
} // for i
} // wrap flag
else if (unwrapflag == 1) {
for (int i = 0; i < bufRange[tid]; i++) {
if ((bufLength[tid] - mpifhStringCountPerThread[tid]) < DUMP_BUF_INCREMENT_SIZE) {
mpifh_buffer_line_per_thread[tid] = (char *) realloc(mpifh_buffer_line_per_thread[tid],(mpifhStringCountPerThread[tid]+DUMP_BUF_CHUNK_SIZE) * sizeof(char));
bufLength[tid] = (mpifhStringCountPerThread[tid]+DUMP_BUF_CHUNK_SIZE) * sizeof(char);
}
for (int j = 0; j < size_one; j++) {
double unwrap_coord;
if (j == 0) {
//offset += sprintf(&sbuf[offset],"%f \n",mybuf[m]);
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),"%f \n",mybuf[bufOffset[tid]+m]);
} else if (j == 1) {
// offset += sprintf(&sbuf[offset],"%s \n",typenames[(int) mybuf[m]]);
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),"%s \n",typenames[(int) mybuf[bufOffset[tid]+m]]);
} else if (j >= 2 && j <= 4) {
unwrap_coord = (mybuf[bufOffset[tid]+m] - 0.5)/UNWRAPEXPAND + 0.5;
//offset += sprintf(&sbuf[offset],vformat[j],unwrap_coord);
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],unwrap_coord);
} else if (j >= 5 ) {
if (vtype[j] == INT)
//offset +=
// sprintf(&sbuf[offset],vformat[j],static_cast<int> (mybuf[m]));
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],static_cast<int> (mybuf[bufOffset[tid]+m]));
else if (vtype[j] == DOUBLE)
// offset += sprintf(&sbuf[offset],vformat[j],mybuf[m]);
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],mybuf[bufOffset[tid]+m]);
else if (vtype[j] == STRING)
// offset +=
// sprintf(&sbuf[offset],vformat[j],typenames[(int) mybuf[m]]);
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],typenames[(int) mybuf[bufOffset[tid]+m]]);
else if (vtype[j] == BIGINT)
// offset +=
// sprintf(&sbuf[offset],vformat[j],static_cast<bigint> (mybuf[m]));
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),vformat[j],static_cast<bigint> (mybuf[bufOffset[tid]+m]));
}
m++;
} // for j
mpifhStringCountPerThread[tid] += sprintf(&(mpifh_buffer_line_per_thread[tid][mpifhStringCountPerThread[tid]]),"\n");
} // for i
} // unwrap flag
} // pragma omp parallel
#pragma omp barrier
mpifhStringCount = 0;
for (i=0;i<nthreads;i++) {
mpifhStringCount += mpifhStringCountPerThread[i];
}
memory->destroy(bufOffset);
memory->destroy(bufRange);
memory->destroy(bufLength);
if (mpifhStringCount > 0) {
if (mpifhStringCount > maxsbuf) {
if (mpifhStringCount > MAXSMALLINT) return -1;
maxsbuf = mpifhStringCount;
memory->grow(sbuf,maxsbuf,"dump:sbuf");
}
sbuf[0] = '\0';
}
for (int i=0;i<nthreads;i++) {
strcat(sbuf,mpifh_buffer_line_per_thread[i]);
free(mpifh_buffer_line_per_thread[i]);
}
memory->destroy(mpifhStringCountPerThread);
free(mpifh_buffer_line_per_thread);
} // else omp
return mpifhStringCount;
}
#endif
/* ---------------------------------------------------------------------- */
void DumpCFGMPIIO::write_data(int n, double *mybuf)
{
(this->*write_choice)(n,mybuf);
}
/* ---------------------------------------------------------------------- */
void DumpCFGMPIIO::write_string(int n, double *mybuf)
{
MPI_Status mpiStatus;
if (performEstimate) {
#if defined(_OPENMP)
int nthreads = omp_get_max_threads();
if (nthreads > 1)
nsme = convert_string_omp(n,mybuf);
else
nsme = convert_string(n,mybuf);
#else
nsme = convert_string(n,mybuf);
#endif
bigint incPrefix = 0;
bigint bigintNsme = (bigint) nsme;
MPI_Scan(&bigintNsme,&incPrefix,1,MPI_LMP_BIGINT,MPI_SUM,world);
sumFileSize = (incPrefix*sizeof(char));
offsetFromHeader = ((incPrefix-bigintNsme)*sizeof(char));
}
else {
MPI_File_write_at_all(mpifh,mpifo+offsetFromHeader,sbuf,nsme,MPI_CHAR,&mpiStatus);
if (flush_flag)
MPI_File_sync(mpifh);
}
}

80
src/MPIIO/dump_cfg_mpiio.h Normal file
View File

@ -0,0 +1,80 @@
/* ----------------------------------------------------------------------
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(cfg/mpiio,DumpCFGMPIIO)
#else
#ifndef LMP_DUMP_CFG_MPIIO_H
#define LMP_DUMP_CFG_MPIIO_H
#include "dump_cfg.h"
namespace LAMMPS_NS {
class DumpCFGMPIIO : public DumpCFG {
public:
DumpCFGMPIIO(class LAMMPS *, int, char **);
virtual ~DumpCFGMPIIO();
protected:
bigint sumFileSize; // size in bytes of the file up through this rank offset from the end of the header data
char *headerBuffer; // buffer for holding header data
MPI_File mpifh;
MPI_Offset mpifo,offsetFromHeader,headerSize, currentFileSize;
int performEstimate; // switch for write_data and write_header methods to use for gathering data and detemining filesize for preallocation vs actually writing the data
char *filecurrent; // name of file for this round (with % and * replaced)
#if defined(_OPENMP)
int convert_string_omp(int, double *); // multithreaded version of convert_string
#endif
virtual void openfile();
virtual void init_style();
virtual void write_header(bigint);
virtual void write();
virtual void write_data(int, double *);
typedef void (DumpCFGMPIIO::*FnPtrData)(int, double *);
FnPtrData write_choice; // ptr to write data functions
void write_string(int, double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Dump cfg arguments must start with 'mass type xs ys zs' or 'mass type xsu ysu zsu'
This is a requirement of the CFG output format. See the dump cfg doc
page for more details.
E: Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu
Self-explanatory.
E: Invalid keyword in dump cfg command
Self-explanatory.
E: Dump cfg requires one snapshot per file
Use the wildcard "*" character in the filename.
*/

1
src/USER-MISC/README
View File

@ -52,4 +52,5 @@ pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13
pair_style lj/sf, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12
pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12
pair_style srp, Tim Sirk, tim.sirk at us.army.mil, 21 Nov 14
pair_style tersoff/table, Luca Ferraro, luca.ferraro@caspur.it, 1 Dec 11

587
src/USER-MISC/fix_srp.cpp Normal file
View File

@ -0,0 +1,587 @@
/* ----------------------------------------------------------------------
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 authors: Timothy Sirk (ARL), Pieter in't Veld (BASF)
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "fix_srp.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
#include "atom_vec.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
// settings
nevery=1;
peratom_freq = 1;
time_integrate = 0;
create_attribute = 0;
comm_border = 2;
// restart settings
restart_global = 1;
restart_peratom = 1;
restart_pbc = 1;
// per-atom array width 2
peratom_flag = 1;
size_peratom_cols = 2;
// initial allocation of atom-based array
// register with Atom class
array = NULL;
grow_arrays(atom->nmax);
// extends pack_exchange()
atom->add_callback(0);
atom->add_callback(1); // restart
atom->add_callback(2);
// zero
for (int i = 0; i < atom->nmax; i++)
for (int m = 0; m < 3; m++)
array[i][m] = 0.0;
// assume setup of fix is needed to insert particles
// yes if reading datafile
// maybe if reading restart
// a datafile cannot contain BPs
// a restart file written during the run has BPs as per atom data
// a restart file written after the run does not have BPs
setup = 1;
}
/* ---------------------------------------------------------------------- */
FixSRP::~FixSRP()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
atom->delete_callback(id,1); // for restart
atom->delete_callback(id,2);
memory->destroy(array);
}
/* ---------------------------------------------------------------------- */
int FixSRP::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
mask |= PRE_EXCHANGE;
mask |= POST_RUN;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSRP::init()
{
if (force->pair_match("hybrid",1) == NULL)
error->all(FLERR,"Cannot use pair srp without pair_style hybrid");
// reserve the highest numbered atom type for bond particles (BPs)
// set bptype if not read from restart
if(!restart_reset) bptype = atom->ntypes;
// check for prescense of extra atom type for bond particle
// if reading a rst file written in the middle of run
// then bond particles are already present
// otherwise need to insert particles in setup
for(int i=0; i < atom->nlocal; i++)
if(atom->type[i] == bptype){
// error if missing extra atom type in either
// data file or rst file without this fix
if(!restart_reset)
error->all(FLERR,"Fix SRP requires an extra atom type");
else
setup = 0;
}
// setup neigh exclusions for diff atom types
// BPs do not interact with other types
// type bptype only interacts with itself
char* arg1[4];
arg1[0] = "exclude";
arg1[1] = "type";
char c0[20];
char c1[20];
for(int z = 1; z < bptype; z++)
{
sprintf(c0, "%d", z);
arg1[2] = c0;
sprintf(c1, "%d", bptype);
arg1[3] = c1;
neighbor->modify_params(4, arg1);
}
}
/* ----------------------------------------------------------------------
insert bond particles
------------------------------------------------------------------------- */
void FixSRP::setup_pre_force(int zz)
{
if(!setup) return;
double rsqold = 0.0;
double delx, dely, delz, rmax, rsq, rsqmax;
double **x = atom->x;
bigint nall = atom->nlocal + atom->nghost;
double xold[nall][3];
// make a copy of coordinates and tags
// create_atom overwrites ghost atoms
for(int i = 0; i < nall; i++){
xold[i][0] = x[i][0];
xold[i][1] = x[i][1];
xold[i][2] = x[i][2];
}
int *tag = atom->tag;
int tagold[nall];
for(int i = 0; i < nall; i++){
tagold[i]=tag[i];
}
int nlocal = atom->nlocal;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int nadd = 0;
int i,j;
for (int n = 0; n < nbondlist; n++) {
// consider only the user defined bond type
if(bondlist[n][2] != btype) continue;
i = bondlist[n][0];
j = bondlist[n][1];
// position of bond i
xone[0] = (xold[i][0] + xold[j][0])*0.5;
xone[1] = (xold[i][1] + xold[j][1])*0.5;
xone[2] = (xold[i][2] + xold[j][2])*0.5;
// record longest bond for ghostcut
delx = xold[j][0] - xold[i][0];
dely = xold[j][1] - xold[i][1];
delz = xold[j][2] - xold[i][2];
rsq = delx*delx + dely*dely + delz*delz;
if(rsq > rsqold) rsqold = rsq;
// make one particle for each bond
// i is local
// if newton bond, always make particle
// if j is local, always make particle
// if j is ghost, decide from tag
if( force->newton_bond || j < nlocal || tagold[i] > tagold[j] ){
atom->natoms += 1;
atom->avec->create_atom(bptype,xone);
// pack tag i/j into buffer for comm
array[atom->nlocal-1][0] = (double)tagold[i];
array[atom->nlocal-1][1] = (double)tagold[j];
nadd++;
}
}
// new total # of atoms
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
// assign tags for new atoms, update map
if (atom->tag_enable) atom->tag_extend();
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
char str[128];
int Nadd = 0;
MPI_Allreduce(&nadd,&Nadd,1,MPI_INT,MPI_SUM,world);
if(comm->me == 0){
sprintf(str, "Inserted %d bond particles.", Nadd);
error->message(FLERR,str);
}
// check ghost comm distances
// warn and change if shorter from estimate
// ghost atoms must be present for bonds on edge of neighbor cutoff
// extend cutghost slightly more than half of the longest bond
MPI_Allreduce(&rsqold,&rsqmax,1,MPI_DOUBLE,MPI_MAX,world);
rmax = sqrt(rsqmax);
double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax;
// find smallest cutghost
double cutghostmin = comm->cutghost[0];
if (cutghostmin > comm->cutghost[1])
cutghostmin = comm->cutghost[1];
if (cutghostmin > comm->cutghost[2])
cutghostmin = comm->cutghost[2];
// reset cutghost if needed
if (cutneighmax_srp > cutghostmin){
if(comm->me == 0){
sprintf(str, "Extending ghost comm cutoff. New %f, old %f.", cutneighmax_srp, cutghostmin);
error->message(FLERR,str);
}
// cutghost updated by comm->setup
comm->cutghostuser = cutneighmax_srp;
}
// put new particles in the box before exchange
// move owned to new procs
// get ghosts
// build neigh lists again
domain->pbc();
comm->setup();
comm->exchange();
if (atom->sortfreq > 0) atom->sort();
comm->borders();
neighbor->build();
neighbor->ncalls = 0;
// new atom counts
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
// zero all forces
for(int i = 0; i < nall; i++)
for(int n = 0; n < 3; n++)
atom->f[i][n] = 0.0;
// do not include bond particles in thermo output
// remove them from all groups
for(int i=0; i< nlocal; i++)
if(atom->type[i] == bptype)
atom->mask[i] = 0;
}
/* ----------------------------------------------------------------------
set position of bond particles
------------------------------------------------------------------------- */
void FixSRP::pre_exchange()
{
// update ghosts
comm->forward_comm();
// reassign bond particle coordinates to midpoint of bonds
// only needed before neigh rebuild
double **x=atom->x;
int i,j;
int nlocal = atom->nlocal;
for(int ii = 0; ii < nlocal; ii++){
if(atom->type[ii] != bptype) continue;
i = atom->map((int)array[ii][0]);
if(i < 0) error->all(FLERR,"Fix SRP failed to map atom");
i = domain->closest_image(ii,i);
j = atom->map((int)array[ii][1]);
if(j < 0) error->all(FLERR,"Fix SRP failed to map atom");
j = domain->closest_image(ii,j);
// position of bond particle ii
atom->x[ii][0] = (x[i][0] + x[j][0])*0.5;
atom->x[ii][1] = (x[i][1] + x[j][1])*0.5;
atom->x[ii][2] = (x[i][2] + x[j][2])*0.5;
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixSRP::memory_usage()
{
double bytes = atom->nmax*2 * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array
------------------------------------------------------------------------- */
void FixSRP::grow_arrays(int nmax)
{
memory->grow(array,nmax,2,"fix_srp:array");
array_atom = array;
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
called when move to new proc
------------------------------------------------------------------------- */
void FixSRP::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < 2; m++)
array[j][m] = array[i][m];
}
/* ----------------------------------------------------------------------
initialize one atom's array values
called when atom is created
------------------------------------------------------------------------- */
void FixSRP::set_arrays(int i)
{
array[i][0] = -1;
array[i][1] = -1;
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixSRP::pack_exchange(int i, double *buf)
{
for (int m = 0; m < 2; m++) buf[m] = array[i][m];
return 2;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixSRP::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < 2; m++) array[nlocal][m] = buf[m];
return 2;
}
/* ----------------------------------------------------------------------
pack values for border communication at re-neighboring
------------------------------------------------------------------------- */
int FixSRP::pack_border(int n, int *list, double *buf)
{
// pack buf for border com
int i,j;
int m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = array[j][0];
buf[m++] = array[j][1];
}
return m;
}
/* ----------------------------------------------------------------------
unpack values for border communication at re-neighboring
------------------------------------------------------------------------- */
int FixSRP::unpack_border(int n, int first, double *buf)
{
// unpack buf into array
int i,last;
int m = 0;
last = first + n;
for (i = first; i < last; i++){
array[i][0] = static_cast<int> (buf[m++]);
array[i][1] = static_cast<int> (buf[m++]);
}
return m;
}
/* ----------------------------------------------------------------------
remove particles after run
------------------------------------------------------------------------- */
void FixSRP::post_run()
{
// all bond particles are removed after each run
// useful for write_data and write_restart commands
// since it occurs between runs
bigint natoms_previous = atom->natoms;
int nlocal = atom->nlocal;
int* dlist;
memory->create(dlist,nlocal,"fix_srp:dlist");
for (int i = 0; i < nlocal; i++){
if(atom->type[i] == bptype)
dlist[i] = 1;
else
dlist[i] = 0;
}
// delete local atoms flagged in dlist
// reset nlocal
AtomVec *avec = atom->avec;
int i = 0;
while (i < nlocal) {
if (dlist[i]) {
avec->copy(nlocal-1,i,1);
dlist[i] = dlist[nlocal-1];
nlocal--;
} else i++;
}
atom->nlocal = nlocal;
memory->destroy(dlist);
// reset atom->natoms
// reset atom->map if it exists
// set nghost to 0 so old ghosts won't be mapped
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// print before and after atom count
bigint ndelete = natoms_previous - atom->natoms;
if (comm->me == 0) {
if (screen) fprintf(screen,"Deleted " BIGINT_FORMAT
" atoms, new total = " BIGINT_FORMAT "\n",
ndelete,atom->natoms);
if (logfile) fprintf(logfile,"Deleted " BIGINT_FORMAT
" atoms, new total = " BIGINT_FORMAT "\n",
ndelete,atom->natoms);
}
// verlet calls box_too_small_check() in post_run
// this check maps all bond partners
// therefore need ghosts
domain->pbc();
comm->setup();
comm->exchange();
if (atom->sortfreq > 0) atom->sort();
comm->borders();
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixSRP::pack_restart(int i, double *buf)
{
int m = 0;
buf[m++] = 3;
buf[m++] = array[i][0];
buf[m++] = array[i][1];
return m;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixSRP::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
int m = 0;
for (int i = 0; i < nth; i++){
m += static_cast<int> (extra[nlocal][m]);
}
m++;
array[nlocal][0] = extra[nlocal][m++];
array[nlocal][1] = extra[nlocal][m++];
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixSRP::maxsize_restart()
{
return 3;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixSRP::size_restart(int nlocal)
{
return 3;
}
/* ----------------------------------------------------------------------
pack global state of Fix
------------------------------------------------------------------------- */
void FixSRP::write_restart(FILE *fp)
{
int n = 0;
double list[3];
list[n++] = comm->cutghostuser;
list[n++] = btype;
list[n++] = bptype;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixSRP::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
comm->cutghostuser = static_cast<double> (list[n++]);
btype = static_cast<int> (list[n++]);
bptype = static_cast<int> (list[n++]);
}
/* ----------------------------------------------------------------------
interface with other classes
pair srp sets the bond type for this fix
------------------------------------------------------------------------- */
int FixSRP::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"btype") == 0) {
btype = atoi(arg[1]);
return 2;
}
return 0;
}

72
src/USER-MISC/fix_srp.h Normal file
View File

@ -0,0 +1,72 @@
/* ----------------------------------------------------------------------
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 FIX_CLASS
FixStyle(SRP,FixSRP)
#else
#ifndef LMP_FIX_SRP_H
#define LMP_FIX_SRP_H
#include "stdio.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixSRP : public Fix {
public:
FixSRP(class LAMMPS *, int, char **);
~FixSRP();
int setmask();
void init();
void pre_exchange();
void setup_pre_force(int);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int, int);
void set_arrays(int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
int pack_border(int, int *, double *);
int unpack_border(int, int, double *);
void post_run();
int pack_restart(int, double*);
void unpack_restart(int, int);
int maxsize_restart();
int size_restart(int);
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
double **array;
private:
double xone[3];
int btype;
int bptype;
int setup;
};
}
#endif
#endif
/* ERROR/WARNING messages:
*/

701
src/USER-MISC/pair_srp.cpp Normal file
View File

@ -0,0 +1,701 @@
/* ----------------------------------------------------------------------
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 authors: Timothy Sirk (ARL), Pieter in't Veld (BASF)
This pair style srp command calculates a segmental repulsive force
between bonds. This is useful for preventing the crossing of bonds if
soft non-bonded potentials are used, such as DPD polymer chains.
See the doc page for pair_style srp command for usage instructions.
There is an example script for this package in examples/USER/srp.
Please contact Timothy Sirk for questions (tim.sirk@us.army.mil).
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "pair_srp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "fix_srp.h"
#include "thermo.h"
#include "output.h"
#include "string.h"
#include "citeme.h"
using namespace LAMMPS_NS;
#define SMALL 1.0e-3
#define BIG 1e10
#define ONETWOBIT 0x40000000
static const char cite_srp[] =
"@Article{Sirk2012\n"
" author = {T. Sirk and Y. Sliozberg and J. Brennan and M. Lisal and J. Andzelm},\n"
" title = {An enhanced entangled polymer model for dissipative particle dynamics},\n"
" journal = {J.~Chem.~Phys.},\n"
" year = 2012,\n"
" volume = 136,\n"
" pages = {134903}\n"
"}\n\n";
/* ----------------------------------------------------------------------
set size of pair comms in constructor
---------------------------------------------------------------------- */
PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp)
{
writedata = 1;
if (lmp->citeme) lmp->citeme->add(cite_srp);
nextra = 1;
segment = NULL;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSRP::allocate()
{
allocated = 1;
// particles of bptype inserted by fix srp
// bptype is the highest numbered atom type
int n = bptype;
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(cut, n + 1, n + 1, "pair:cut");
memory->create(a0, n + 1, n + 1, "pair:a0");
// setflag for atom types
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
maxcount = 0;
}
/* ----------------------------------------------------------------------
free
------------------------------------------------------------------------- */
PairSRP::~PairSRP()
{
if (allocated)
{
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(a0);
memory->destroy(segment);
}
// check nfix in case all fixes have already been deleted
if (modify->nfix) modify->delete_fix("mysrp");
}
/* ----------------------------------------------------------------------
compute bond-bond repulsions
------------------------------------------------------------------------- */
void PairSRP::compute(int eflag, int vflag)
{
// setup energy and virial
if (eflag || vflag)
ev_setup(eflag, vflag);
else
evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int i0, i1, j0, j1;
int i,j,ii,jj,inum,jnum;
double dijsq, dij;
int *ilist,*jlist,*numneigh,**firstneigh;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double dx,dy,dz,ti,tj;
double wd, lever0, lever1, evdwl, fpair;
double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1;
double fx, fy, fz;
// mapping global to local for atoms inside bond particles
// exclude 1-2 neighs if requested
if (neighbor->ago == 0){
remapBonds(nall);
if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh);
}
// this pair style only used with hybrid
// so each atom i is type bptype
// each neigh j is type bptype due to exclusions
// using midpoint distance option
if(midpoint){
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jnum = numneigh[i];
// two atoms inside bond particle
i0 = segment[i][0];
j0 = segment[i][1];
for (jj = 0; jj < jnum; jj++) {
jlist = firstneigh[i];
j = jlist[jj];
// enforce 1-2 exclusions
if( (sbmask(j) & exclude) )
continue;
j &= NEIGHMASK;
//atoms inside bond particle
i1 = segment[j][0];
j1 = segment[j][1];
// midpt dist bond 0 and 1
dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]);
dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]);
dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]);
dijsq = dx*dx + dy*dy + dz*dz;
if (dijsq < cutsq[bptype][bptype]){
dij = pow(dijsq, 0.5);
if (dij < SMALL)
dij = SMALL; // prevent explosions
wd = 1.0 - dij / cut[bptype][bptype];
fpair = 0.5 * a0[bptype][bptype] * wd / dij;
// force for bond 0, beads 0,1
//force between bonds
fx = fpair * dx;
fy = fpair * dy;
fz = fpair * dz;
f[i0][0] += fx; //keep force sign for bond 0
f[i0][1] += fy;
f[i0][2] += fz;
f[j0][0] += fx;
f[j0][1] += fy;
f[j0][2] += fz;
f[i1][0] -= fx; //flip force sign for bond 1
f[i1][1] -= fy;
f[i1][2] -= fz;
f[j1][0] -= fx;
f[j1][1] -= fy;
f[j1][2] -= fz;
// ************************************************* //
if (eflag){
evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd;
}
if (evflag){
ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz);
ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz);
}
if (vflag_fdotr) virial_fdotr_compute();
}
}
}
}
else{
// using min distance option
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jnum = numneigh[i];
i0 = segment[i][0];
j0 = segment[i][1];
for (jj = 0; jj < jnum; jj++) {
jlist = firstneigh[i];
j = jlist[jj];
// enforce 1-2 exclusions
if( (sbmask(j) & exclude) )
continue;
j &= NEIGHMASK;
i1 = segment[j][0];
j1 = segment[j][1];
getMinDist(x, dx, dy, dz, ti, tj, i0, j0, i1, j1);
dijsq = dx*dx + dy*dy + dz*dz;
if (dijsq < cutsq[bptype][bptype]){
dij = pow(dijsq, 0.5);
if (dij < SMALL)
dij = SMALL; // prevent explosions
wd = 1.0 - dij / cut[bptype][bptype];
fpair = a0[bptype][bptype] * wd / dij;
// force for bond 0, beads 0,1
lever0 = 0.5 + ti; // assign force according to lever rule
lever1 = 0.5 + tj; // assign force according to lever rule
//force between bonds
fx = fpair * dx;
fy = fpair * dy;
fz = fpair * dz;
//decompose onto atoms
fxlever0 = fx * lever0;
fylever0 = fy * lever0;
fzlever0 = fz * lever0;
fxlever1 = fx * lever1;
fylever1 = fy * lever1;
fzlever1 = fz * lever1;
f[i0][0] += fxlever0; //keep force sign for bond 0
f[i0][1] += fylever0;
f[i0][2] += fzlever0;
f[j0][0] += (fx - fxlever0);
f[j0][1] += (fy - fylever0);
f[j0][2] += (fz - fzlever0);
f[i1][0] -= fxlever1; //flip force sign for bond 1
f[i1][1] -= fylever1;
f[i1][2] -= fzlever1;
f[j1][0] -= (fx - fxlever1);
f[j1][1] -= (fy - fylever1);
f[j1][2] -= (fz - fzlever1);
// ************************************************* //
if (eflag){
evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd;
}
if (evflag){
ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz);
ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz);
}
if (vflag_fdotr) virial_fdotr_compute();
}
}
}
}
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSRP::settings(int narg, char **arg)
{
if (narg < 3 || narg > 5)
error->all(FLERR,"Illegal pair_style command");
cut_global = force->numeric(FLERR,arg[0]);
btype = force->inumeric(FLERR,arg[1]);
if (btype > atom->nbondtypes) error->all(FLERR,"Illegal pair_style command");
// settings
midpoint = 0;
min = 0;
if (strcmp(arg[2],"min") == 0) min = 1;
else if (strcmp(arg[2],"mid") == 0) midpoint = 1;
else
error->all(FLERR,"Illegal pair_style command");
int iarg = 3;
// default exclude 1-2
// scaling for 1-2, etc not supported
exclude = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"exclude") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp exclude command");
if (strcmp(arg[iarg+1],"yes") == 0)
exclude = 1;
if (strcmp(arg[iarg+1],"no") == 0)
exclude = 0;
iarg += 2;
} else error->all(FLERR,"Illegal pair srp command");
}
// highest atom type is saved for bond particles (BPs)
bptype = atom->ntypes;
// reset cutoffs if explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= bptype; i++)
for (j = i+1; j <= bptype; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
set coeffs
------------------------------------------------------------------------- */
void PairSRP::coeff(int narg, char **arg)
{
if (narg < 3 || narg > 4)
error->all(FLERR,"PairSRP: Incorrect args for pair coeff");
if (!allocated) allocate();
// set ij bond-bond cutoffs
int ilo, ihi, jlo, jhi;
force->bounds(arg[0], bptype, ilo, ihi);
force->bounds(arg[1], bptype, jlo, jhi);
double a0_one = force->numeric(FLERR,arg[2]);
double cut_one = cut_global;
if (narg == 4) cut_one = force->numeric(FLERR,arg[3]);
int count = 0;
for (int i = ilo; i <= ihi; i++)
{
for (int j = MAX(jlo,i); j <= jhi; j++)
{
a0[i][j] = a0_one;
cut[i][j] = cut_one;
cutsq[i][j] = cut_one * cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->warning(FLERR,"PairSRP: No pair coefficients were set");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairSRP::init_style()
{
if (!force->newton_pair)
error->all(FLERR,"PairSRP: Pair srp requires newton pair on");
// need fix srp
// invoke here instead of constructor
// to make restart possible
char **newarg = new char*[3];
newarg[0] = (char *) "pair_srp";
newarg[1] = (char *) "all";
newarg[2] = (char *) "SRP";
modify->add_fix(3,newarg);
f_srp = (FixSRP *) modify->fix[modify->nfix-1];
delete [] newarg;
// set bond type in fix srp
// bonds of this type will be represented by bond particles (BPs)
// btype = bond type
// bptype = bond particle type
char c0[20];
char* arg0[2];
sprintf(c0, "%d", btype);
arg0[0] = (char *) "btype";
arg0[1] = c0;
f_srp->modify_params(2, arg0);
// BPs do not contribute to energy or virial
// BPs do not belong to group all
// but thermo normalization is by nall
// turn off normalization
int me;
MPI_Comm_rank(world,&me);
char *arg1[2];
arg1[0] = (char *) "norm";
arg1[1] = (char *) "no";
output->thermo->modify_params(2, arg1);
if (me == 0)
error->message(FLERR,"Thermo normalization turned off by pair srp");
neighbor->request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSRP::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"PairSRP: All pair coeffs are not set");
cut[j][i] = cut[i][j];
a0[j][i] = a0[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
find min distance for bonds i0/j0 and i1/j1
------------------------------------------------------------------------- */
inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1)
{
// move these outside the loop
double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj;
double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1;
// compute midpt dist from 1st atom, 1st bond
diffx0 = x[j0][0] - x[i0][0]; // x,y,z from bond 0
diffy0 = x[j0][1] - x[i0][1];
diffz0 = x[j0][2] - x[i0][2];
// compute midpt dist from 1st atom, 2nd bond
diffx1 = x[j1][0] - x[i1][0];
diffy1 = x[j1][1] - x[i1][1];
diffz1 = x[j1][2] - x[i1][2];
// midpoint distance
dPx = 0.5*(diffx0-diffx1) + x[i0][0]-x[i1][0];
dPy = 0.5*(diffy0-diffy1) + x[i0][1]-x[i1][1];
dPz = 0.5*(diffz0-diffz1) + x[i0][2]-x[i1][2];
// Ri^2 Rj^2
RiRi = diffx0*diffx0 + diffy0*diffy0 + diffz0*diffz0;
RiRj = diffx0*diffx1 + diffy0*diffy1 + diffz0*diffz1;
RjRj = diffx1*diffx1 + diffy1*diffy1 + diffz1*diffz1;
denom = RiRj*RiRj - RiRi*RjRj;
// handle case of parallel lines
// reduce to midpt distance
if (fabs(denom) < SMALL){
if(denom < 0) denom = -BIG;
else denom = BIG;
}
// calc ti
termx0 = RiRj*diffx1 - RjRj*diffx0;
termy0 = RiRj*diffy1 - RjRj*diffy0;
termz0 = RiRj*diffz1 - RjRj*diffz0;
num0 = dPx*termx0 + dPy*termy0 + dPz*termz0;
ti = num0 / denom;
if (ti > 0.5) ti = 0.5;
if (ti < -0.5) ti = -0.5;
// calc tj
termx1 = RiRj*diffx0 - RiRi*diffx1;
termy1 = RiRj*diffy0 - RiRi*diffy1;
termz1 = RiRj*diffz0 - RiRi*diffz1;
num1 = dPx*termx1 + dPy*termy1 + dPz*termz1;
tj = -num1/ denom;
if (tj > 0.5) tj = 0.5;
if (tj < -0.5) tj = -0.5;
// min dist
dx = dPx - ti*diffx0 + tj*diffx1;
dy = dPy - ti*diffy0 + tj*diffy1;
dz = dPz - ti*diffz0 + tj*diffz1;
}
/* --------------------------------------------------------
map global id of atoms in stored by each bond particle
------------------------------------------------------- */
inline void PairSRP::remapBonds(int &nall)
{
if(nall > maxcount){
memory->grow(segment, nall, 2, "pair:segment");
maxcount = nall;
}
// loop over all bond particles
// each bond paricle holds two bond atoms
// map global ids of bond atoms to local ids
// might not be able to map both bond atoms of j, if j is outside neighcut
// these are not on neighlist, so are not used
int tmp;
srp = f_srp->array_atom;
for (int i = 0; i < nall; i++) {
if(atom->type[i] == bptype){
// tmp is local id
// tmp == -1 is ok
tmp = atom->map((int)srp[i][0]);
segment[i][0] = domain->closest_image(i,tmp);
// repeat with other id
tmp = atom->map((int)srp[i][1]);
segment[i][1] = domain->closest_image(i,tmp);
}
}
}
/* --------------------------------------------------------
add exclusions for 1-2 neighs, if requested
more complex exclusions or scaling probably not needed
------------------------------------------------------- */
inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh)
{
int i0, i1, j0, j1;
int i,j,ii,jj,jnum;
// encode neighs with exclusions
// only need 1-2 info for normal uses of srp
// add 1-3, etc later if ever needed
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jnum = numneigh[i];
// two atoms inside bond particle
i0 = segment[i][0];
j0 = segment[i][1];
for (jj = 0; jj < jnum; jj++) {
jlist = firstneigh[i];
j = jlist[jj];
j &= NEIGHMASK;
//two atoms inside bond particle
i1 = segment[j][0];
j1 = segment[j][1];
// check for a 1-2 neigh
if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){
j |= ONETWOBIT;
jlist[jj] = j;
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairSRP::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d %g\n",i,a0[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairSRP::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]);
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSRP::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&a0[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSRP::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
printf(" i %d j %d \n",i,j);
fread(&a0[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSRP::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&bptype,sizeof(int),1,fp);
fwrite(&btype,sizeof(int),1,fp);
fwrite(&min,sizeof(int),1,fp);
fwrite(&midpoint,sizeof(int),1,fp);
fwrite(&exclude,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSRP::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&bptype,sizeof(int),1,fp);
fread(&btype,sizeof(int),1,fp);
fread(&min,sizeof(int),1,fp);
fread(&midpoint,sizeof(int),1,fp);
fread(&exclude,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
}

68
src/USER-MISC/pair_srp.h Normal file
View File

@ -0,0 +1,68 @@
/* ----------------------------------------------------------------------
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 PAIR_CLASS
PairStyle(srp,PairSRP)
#else
#ifndef LMP_PAIR_SRP_H
#define LMP_PAIR_SRP_H
#include "pair.h"
namespace LAMMPS_NS {
class PairSRP : public Pair {
public:
PairSRP(class LAMMPS *);
virtual ~PairSRP();
virtual void compute(int, int);
virtual void settings(int, char **);
virtual void coeff(int, char **);
void init_style();
double init_one(int, int);
virtual void write_data(FILE *);
virtual void write_data_all(FILE *);
virtual void write_restart(FILE *);
virtual void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
private:
inline void onetwoexclude(int* &, int &, int* &, int* &, int** &);
inline void remapBonds(int &);
void allocate();
void getMinDist(double** &, double &, double &, double &, double &,
double &, int &, int &, int &, int &);
bool min, midpoint;
double **cut;
double **a0;
double **srp;
double cut_global;
int bptype;
int btype;
class Fix *f_srp;
int exclude,maxcount;
int **segment;
};
}
#endif
#endif
/* ERROR/WARNING messages:
*/