From 13c55490099890791a2bed3aff92ced293a53da8 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 7 Apr 2016 21:12:44 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14809 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/bond.cpp | 71 ++++++++++++++++++++++++++++++++++++++++++++++ src/bond.h | 2 ++ src/input.cpp | 12 ++++++++ src/input.h | 1 + src/pair.cpp | 8 +++--- src/pair_table.cpp | 27 +++++++++++++----- 6 files changed, 110 insertions(+), 11 deletions(-) diff --git a/src/bond.cpp b/src/bond.cpp index 443816418f..5c2622281d 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -16,6 +16,7 @@ #include "atom.h" #include "comm.h" #include "force.h" +#include "neighbor.h" #include "suffix.h" #include "atom_masks.h" #include "memory.h" @@ -23,6 +24,8 @@ using namespace LAMMPS_NS; +enum{NONE,LINEAR,SPLINE}; + /* ----------------------------------------------------------------------- set bond contribution to Vdwl energy to 0.0 a particular bond style can override this @@ -212,6 +215,74 @@ void Bond::ev_tally(int i, int j, int nlocal, int newton_bond, } } +/* ---------------------------------------------------------------------- + write a table of bond potential energy/force vs distance to a file +------------------------------------------------------------------------- */ + +void Bond::write_file(int narg, char **arg) +{ + if (narg != 6 && narg !=8) error->all(FLERR,"Illegal bond_write command"); + + // parse optional arguments + + int itype = 0; + int jtype = 0; + if (narg == 8) { + itype = force->inumeric(FLERR,arg[6]); + jtype = force->inumeric(FLERR,arg[7]); + if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes) + error->all(FLERR,"Invalid atom types in bond_write command"); + } + + int btype = force->inumeric(FLERR,arg[0]); + int n = force->inumeric(FLERR,arg[1]); + double inner = force->numeric(FLERR,arg[2]); + double outer = force->numeric(FLERR,arg[3]); + if (inner <= 0.0 || inner >= outer) + error->all(FLERR,"Invalid rlo/rhi values in bond_write command"); + + + double r0 = equilibrium_distance(btype); + + // open file in append mode + // print header in format used by bond_style table + + int me; + MPI_Comm_rank(world,&me); + FILE *fp; + if (me == 0) { + fp = fopen(arg[4],"a"); + if (fp == NULL) error->one(FLERR,"Cannot open bond_write file"); + } + + // initialize potentials before evaluating bond potential + // insures all bond coeffs are set and force constants + // also initialize neighbor so that neighbor requests are processed + // NOTE: might be safest to just do lmp->init() + + force->init(); + neighbor->init(); + + if (me == 0) { + double r,e,f; + + // evaluate energy and force at each of N distances + // note that Bond::single() takes r**2 and returns f/r. + + fprintf(fp,"# Bond potential %s for bond type %d: i,r,energy,force\n", + force->bond_style,btype); + fprintf(fp,"\n%s\nN %d EQ %.15g\n\n",arg[5],n,r0); + + const double dr = (outer-inner) / static_cast(n-1); + for (int i = 0; i < n; i++) { + r = inner + dr * static_cast(i); + e = single(btype,r*r,itype,jtype,f); + fprintf(fp,"%d %.15g %.15g %.15g\n",i+1,r,e,f*r); + } + fclose(fp); + } +} + /* ---------------------------------------------------------------------- */ double Bond::memory_usage() diff --git a/src/bond.h b/src/bond.h index 6841e9362d..d455cda204 100644 --- a/src/bond.h +++ b/src/bond.h @@ -54,6 +54,8 @@ class Bond : protected Pointers { virtual unsigned int data_mask() {return datamask;} virtual unsigned int data_mask_ext() {return datamask_ext;} + void write_file(int, char**); + protected: int suffix_flag; // suffix compatibility flag diff --git a/src/input.cpp b/src/input.cpp index 3ceb906911..bcb64effe9 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -644,6 +644,7 @@ int Input::execute_command() else if (!strcmp(command,"atom_style")) atom_style(); else if (!strcmp(command,"bond_coeff")) bond_coeff(); else if (!strcmp(command,"bond_style")) bond_style(); + else if (!strcmp(command,"bond_write")) bond_write(); else if (!strcmp(command,"boundary")) boundary(); else if (!strcmp(command,"box")) box(); else if (!strcmp(command,"comm_modify")) comm_modify(); @@ -1281,6 +1282,17 @@ void Input::bond_style() /* ---------------------------------------------------------------------- */ +void Input::bond_write() +{ + if (atom->avec->bonds_allow == 0) + error->all(FLERR,"Bond_write command when no bonds allowed"); + if (force->bond == NULL) + error->all(FLERR,"Bond_write command before bond_style is defined"); + else force->bond->write_file(narg,arg); +} + +/* ---------------------------------------------------------------------- */ + void Input::boundary() { if (domain->box_exist) diff --git a/src/input.h b/src/input.h index 9bdedabe14..ccda3f49b1 100644 --- a/src/input.h +++ b/src/input.h @@ -86,6 +86,7 @@ class Input : protected Pointers { void atom_style(); void bond_coeff(); void bond_style(); + void bond_write(); void boundary(); void box(); void comm_modify(); diff --git a/src/pair.cpp b/src/pair.cpp index 89c5bea380..e3c4a1aa36 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1574,9 +1574,9 @@ void Pair::write_file(int narg, char **arg) fprintf(fp,"# Pair potential %s for atom types %d %d: i,r,energy,force\n", force->pair_style,itype,jtype); if (style == RLINEAR) - fprintf(fp,"\n%s\nN %d R %g %g\n\n",arg[7],n,inner,outer); + fprintf(fp,"\n%s\nN %d R %.15g %.15g\n\n",arg[7],n,inner,outer); if (style == RSQ) - fprintf(fp,"\n%s\nN %d RSQ %g %g\n\n",arg[7],n,inner,outer); + fprintf(fp,"\n%s\nN %d RSQ %.15g %.15g\n\n",arg[7],n,inner,outer); } // initialize potentials before evaluating pair potential @@ -1618,7 +1618,7 @@ void Pair::write_file(int narg, char **arg) init_bitmap(inner,outer,n,masklo,maskhi,nmask,nshiftbits); int ntable = 1 << n; if (me == 0) - fprintf(fp,"\n%s\nN %d BITMAP %g %g\n\n",arg[7],ntable,inner,outer); + fprintf(fp,"\n%s\nN %d BITMAP %.15g %.15g\n\n",arg[7],ntable,inner,outer); n = ntable; } @@ -1647,7 +1647,7 @@ void Pair::write_file(int narg, char **arg) e = single(0,1,itype,jtype,rsq,1.0,1.0,f); f *= r; } else e = f = 0.0; - if (me == 0) fprintf(fp,"%d %g %g %g\n",i+1,r,e,f); + if (me == 0) fprintf(fp,"%d %.15g %.15g %.15g\n",i+1,r,e,f); } // restore original vecs that were swapped in for diff --git a/src/pair_table.cpp b/src/pair_table.cpp index 1013bd3677..71b4d49b5f 100644 --- a/src/pair_table.cpp +++ b/src/pair_table.cpp @@ -395,13 +395,16 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) union_int_float_t rsq_lookup; int rerror = 0; + int cerror = 0; fgets(line,MAXLINE,fp); for (int i = 0; i < tb->ninput; i++) { - fgets(line,MAXLINE,fp); - sscanf(line,"%d %lg %lg %lg",&itmp,&rfile,&tb->efile[i],&tb->ffile[i]); - rnew = rfile; + if (NULL == fgets(line,MAXLINE,fp)) + error->one(FLERR,"Premature end of file in pair table"); + if (4 != sscanf(line,"%d %lg %lg %lg", + &itmp,&rfile,&tb->efile[i],&tb->ffile[i])) ++cerror; + rnew = rfile; if (tb->rflag == RLINEAR) rnew = tb->rlo + (tb->rhi - tb->rlo)*i/(tb->ninput-1); else if (tb->rflag == RSQ) { @@ -419,6 +422,7 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) } if (tb->rflag && fabs(rnew-rfile)/rfile > EPSILONR) rerror++; + tb->rfile[i] = rnew; } @@ -452,8 +456,8 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) if (ferror) { char str[128]; - sprintf(str,"%d force values in table are inconsistent with -dE/dr; " - "should only be mistakenly flagged at inflection points",ferror); + sprintf(str,"%d of %d force values in table are inconsistent with -dE/dr.\n" + " Should only be flagged at inflection points",ferror,tb->ninput); error->warning(FLERR,str); } @@ -461,8 +465,17 @@ void PairTable::read_table(Table *tb, char *file, char *keyword) if (rerror) { char str[128]; - sprintf(str,"%d distance values in table differ signifcantly " - "from re-computed values",rerror); + sprintf(str,"%d of %d distance values in table with relative error\n" + " over %g to re-computed values",rerror,tb->ninput,EPSILONR); + error->warning(FLERR,str); + } + + // warn if data was read incompletely, e.g. columns were missing + + if (cerror) { + char str[128]; + sprintf(str,"%d of %d lines in table were incomplete\n" + " or could not be parsed completely",cerror,tb->ninput); error->warning(FLERR,str); } }