From 658d2c76695d708c6dece1f8f9ef6774d60f69b6 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 10 Nov 2006 21:29:18 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@136 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/dump_custom.cpp | 229 ++++++++++++++++++++++++++++++++------------ src/dump_custom.h | 38 ++++++-- src/thermo.cpp | 20 ++-- src/thermo.h | 2 +- 4 files changed, 207 insertions(+), 82 deletions(-) diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 6bea3257f7..c3eb621fa9 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "dump_custom.h" #include "atom.h" +#include "force.h" #include "domain.h" #include "region.h" #include "group.h" @@ -30,7 +31,8 @@ // customize by adding to 1st enum enum{TAG,MOL,TYPE,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ,VX,VY,VZ,FX,FY,FZ, - Q,MUX,MUY,MUZ,TQX,TQY,TQZ,CENTRO,ENG,SXX,SYY,SZZ,SXY,SXZ,SYZ}; + Q,MUX,MUY,MUZ,TQX,TQY,TQZ, + ETOTAL,KE,EPAIR,CENTRO,SXX,SYY,SZZ,SXY,SXZ,SYZ}; enum{LT,LE,GT,GE,EQ,NEQ}; enum{INT,DOUBLE}; @@ -38,6 +40,8 @@ enum{INT,DOUBLE}; DumpCustom::DumpCustom(int narg, char **arg) : Dump(narg, arg) { + int i; + if (narg == 5) error->all("No dump custom arguments specified"); nevery = atoi(arg[3]); @@ -52,13 +56,16 @@ DumpCustom::DumpCustom(int narg, char **arg) : Dump(narg, arg) thresh_op = NULL; thresh_value = NULL; + fixflag_epair = 0; fixflag_centro = 0; - fixflag_energy = 0; fixflag_stress = 0; + computeflag_etotal = 0; + computeflag_ke = 0; + + // process command args // customize by adding to if statement - int i; for (int iarg = 5; iarg < narg; iarg++) { i = iarg-5; if (strcmp(arg[iarg],"tag") == 0) { @@ -161,14 +168,24 @@ DumpCustom::DumpCustom(int narg, char **arg) : Dump(narg, arg) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqz; vtype[i] = DOUBLE; + } else if (strcmp(arg[iarg],"etotal") == 0) { + pack_choice[i] = &DumpCustom::pack_etotal; + vtype[i] = DOUBLE; + computeflag_etotal = 1; + computeflag_ke = 1; + fixflag_epair = 1; + } else if (strcmp(arg[iarg],"ke") == 0) { + pack_choice[i] = &DumpCustom::pack_ke; + vtype[i] = DOUBLE; + computeflag_ke = 1; + } else if (strcmp(arg[iarg],"epair") == 0) { + pack_choice[i] = &DumpCustom::pack_epair; + vtype[i] = DOUBLE; + fixflag_epair = 1; } else if (strcmp(arg[iarg],"centro") == 0) { pack_choice[i] = &DumpCustom::pack_centro; vtype[i] = DOUBLE; fixflag_centro = 1; - } else if (strcmp(arg[iarg],"eng") == 0) { - pack_choice[i] = &DumpCustom::pack_energy; - vtype[i] = DOUBLE; - fixflag_energy = 1; } else if (strcmp(arg[iarg],"sxx") == 0) { pack_choice[i] = &DumpCustom::pack_sxx; vtype[i] = DOUBLE; @@ -196,10 +213,19 @@ DumpCustom::DumpCustom(int narg, char **arg) : Dump(narg, arg) } else error->all("Invalid keyword in dump custom command"); } + if (fixflag_epair) fix_create("ENERGY"); if (fixflag_centro) fix_create("CENTRO"); - if (fixflag_energy) fix_create("ENERGY"); if (fixflag_stress) fix_create("STRESS"); + maxlocal = 0; + choose = NULL; + dchoose = NULL; + + if (computeflag_etotal) etotal = NULL; + if (computeflag_ke) ke = NULL; + + // setup format strings + vformat = new char*[size_one]; format_default = new char[3*size_one+1]; @@ -211,10 +237,6 @@ DumpCustom::DumpCustom(int narg, char **arg) : Dump(narg, arg) vformat[i] = NULL; } - maxchoose = 0; - choose = NULL; - dchoose = NULL; - // one-time file open if (multifile == 0) openfile(); @@ -226,10 +248,10 @@ DumpCustom::DumpCustom(int narg, char **arg) : Dump(narg, arg) DumpCustom::~DumpCustom() { + delete [] pack_choice; delete [] vtype; for (int i = 0; i < size_one; i++) delete [] vformat[i]; delete [] vformat; - delete [] pack_choice; if (nthresh) { memory->sfree(thresh_array); @@ -239,14 +261,17 @@ DumpCustom::~DumpCustom() // delete fixes for computing per-atom custom quantities if they exist + if (fixflag_epair) fix_delete("ENERGY"); if (fixflag_centro) fix_delete("CENTRO"); - if (fixflag_energy) fix_delete("ENERGY"); if (fixflag_stress) fix_delete("STRESS"); - if (maxchoose) { - delete [] choose; - delete [] dchoose; - } + // delete choose arrays and local compute arrys if they exist + + delete [] choose; + delete [] dchoose; + + if (computeflag_etotal) delete [] etotal; + if (computeflag_ke) delete [] ke; } /* ---------------------------------------------------------------------- */ @@ -285,33 +310,32 @@ void DumpCustom::init() // set fixflag if dump invokes any fixes fixflag = 0; - if (fixflag_centro || fixflag_energy || fixflag_stress) fixflag = 1; + if (fixflag_epair || fixflag_centro || fixflag_stress) fixflag = 1; // find associated fixes that must exist // could have changed locations in fix list since created + if (fixflag_epair) ifix_epair = fix_match("ENERGY"); if (fixflag_centro) ifix_centro = fix_match("CENTRO"); - if (fixflag_energy) ifix_energy = fix_match("ENERGY"); if (fixflag_stress) ifix_stress = fix_match("STRESS"); - // allocate choose arrays if first time + // set computeflag if dump invokes any local computation - if (maxchoose == 0) { - maxchoose = atom->nmax; - choose = new int[maxchoose]; - dchoose = new double[maxchoose]; - } + computeflag = 0; + if (computeflag_etotal || computeflag_ke) computeflag = 1; } /* ---------------------------------------------------------------------- - return # of bytes of allocated memory in buf and choose arrays + return # of bytes of allocated memory in buf and choose and local arrays ------------------------------------------------------------------------- */ int DumpCustom::memory_usage() { int bytes = maxbuf * sizeof(double); - bytes += maxchoose * sizeof(int); - bytes += maxchoose * sizeof(double); + bytes += maxlocal * sizeof(int); + bytes += maxlocal * sizeof(double); + if (computeflag_etotal) bytes += maxlocal * sizeof(double); + if (computeflag_ke) bytes += maxlocal * sizeof(double); return bytes; } @@ -386,8 +410,10 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"tqx") == 0) thresh_array[nthresh] = TQX; else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY; else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ; + else if (strcmp(arg[1],"etotal") == 0) thresh_array[nthresh] = ETOTAL; + else if (strcmp(arg[1],"ke") == 0) thresh_array[nthresh] = KE; + else if (strcmp(arg[1],"epair") == 0) thresh_array[nthresh] = EPAIR; else if (strcmp(arg[1],"centro") == 0) thresh_array[nthresh] = CENTRO; - else if (strcmp(arg[1],"eng") == 0) thresh_array[nthresh] = ENG; else if (strcmp(arg[1],"sxx") == 0) thresh_array[nthresh] = SXX; else if (strcmp(arg[1],"syy") == 0) thresh_array[nthresh] = SYY; else if (strcmp(arg[1],"szz") == 0) thresh_array[nthresh] = SZZ; @@ -408,14 +434,14 @@ int DumpCustom::modify_param(int narg, char **arg) // add new fix if necessary + if (thresh_array[nthresh] == EPAIR && fixflag_epair == 0) { + fix_create("ENERGY"); + fixflag_epair = 1; + } if (thresh_array[nthresh] == CENTRO && fixflag_centro == 0) { fix_create("CENTRO"); fixflag_centro = 1; } - if (thresh_array[nthresh] == ENG && fixflag_energy == 0) { - fix_create("ENERGY"); - fixflag_energy = 1; - } if ((thresh_array[nthresh] == SXX || thresh_array[nthresh] == SYY || thresh_array[nthresh] == SZZ || thresh_array[nthresh] == SXY || thresh_array[nthresh] == SXZ || thresh_array[nthresh] == SYZ) && @@ -443,23 +469,35 @@ int DumpCustom::count() { int i; + // grow choose and local-computation arrays if needed + + int nlocal = atom->nlocal; + if (nlocal > maxlocal) { + delete [] choose; + delete [] dchoose; + if (computeflag_etotal) delete [] etotal; + if (computeflag_ke) delete [] ke; + maxlocal = atom->nmax; + choose = new int[maxlocal]; + dchoose = new double[maxlocal]; + if (computeflag_etotal) etotal = new double[maxlocal]; + if (computeflag_ke) ke = new double[maxlocal]; + } + // invoke any fixes that compute dump quantities if (fixflag) { + if (fixflag_epair) modify->fix[ifix_epair]->dump(); if (fixflag_centro) modify->fix[ifix_centro]->dump(); - if (fixflag_energy) modify->fix[ifix_energy]->dump(); if (fixflag_stress) modify->fix[ifix_stress]->dump(); } - // grow choose and dchoose arrays if needed + // invoke any local computation of dump quantities + // etotal must come after ke and fix_epair since uses their results - int nlocal = atom->nlocal; - if (nlocal > maxchoose) { - delete [] choose; - delete [] dchoose; - maxchoose = atom->nmax; - choose = new int[maxchoose]; - dchoose = new double[maxchoose]; + if (computeflag) { + if (computeflag_ke) compute_ke(); + if (computeflag_etotal) compute_etotal(); } // choose all local atoms for output @@ -631,12 +669,18 @@ int DumpCustom::count() } else if (thresh_array[ithresh] == TQZ) { ptr = &atom->torque[0][2]; nstride = 3; + } else if (thresh_array[ithresh] == ETOTAL) { + ptr = etotal; + nstride = 1; + } else if (thresh_array[ithresh] == KE) { + ptr = ke; + nstride = 1; + } else if (thresh_array[ithresh] == EPAIR) { + ptr = ((FixEnergy *) modify->fix[ifix_epair])->energy; + nstride = 3; } else if (thresh_array[ithresh] == CENTRO) { ptr = ((FixCentro *) modify->fix[ifix_centro])->centro; nstride = 3; - } else if (thresh_array[ithresh] == ENG) { - ptr = ((FixEnergy *) modify->fix[ifix_energy])->energy; - nstride = 3; } else if (thresh_array[ithresh] == SXX) { ptr = &((FixStress *) modify->fix[ifix_stress])->stress[0][0]; nstride = 6; @@ -839,7 +883,7 @@ int DumpCustom::fix_match(char *keyword) } // ---------------------------------------------------------------------- -// one routine for every kind of quantity dump custom can output +// one method for every kind of quantity dump custom can output // the atom quantity is packed into buf starting at n with stride size_one // customize by adding a method // ---------------------------------------------------------------------- @@ -1253,6 +1297,47 @@ void DumpCustom::pack_tqz(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_etotal(int n) +{ + double *epair = ((FixEnergy *) modify->fix[ifix_epair])->energy; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = ke[i] + epair[i]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_ke(int n) +{ + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = ke[i]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_epair(int n) +{ + double *epair = ((FixEnergy *) modify->fix[ifix_epair])->energy; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = epair[i]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_centro(int n) { double *centro = ((FixCentro *) modify->fix[ifix_centro])->centro; @@ -1267,20 +1352,6 @@ void DumpCustom::pack_centro(int n) /* ---------------------------------------------------------------------- */ -void DumpCustom::pack_energy(int n) -{ - double *energy = ((FixEnergy *) modify->fix[ifix_energy])->energy; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (choose[i]) { - buf[n] = energy[i]; - n += size_one; - } -} - -/* ---------------------------------------------------------------------- */ - void DumpCustom::pack_sxx(int n) { double **stress = ((FixStress *) modify->fix[ifix_stress])->stress; @@ -1362,3 +1433,39 @@ void DumpCustom::pack_syz(int n) n += size_one; } } + +// ---------------------------------------------------------------------- +// one method for each locally-computed dump quantity +// ---------------------------------------------------------------------- + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::compute_etotal() +{ + double *epair = ((FixEnergy *) modify->fix[ifix_epair])->energy; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + etotal[i] = ke[i] + epair[i]; +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::compute_ke() +{ + double mvv2e = force->mvv2e; + double **v = atom->v; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int nlocal = atom->nlocal; + + if (atom->mass_require) + for (int i = 0; i < nlocal; i++) + ke[i] = 0.5 * mvv2e * mass[type[i]] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + else + for (int i = 0; i < nlocal; i++) + ke[i] = 0.5 * mvv2e * rmass[i] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); +} diff --git a/src/dump_custom.h b/src/dump_custom.h index c58f965338..fe836d5c69 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -34,20 +34,31 @@ class DumpCustom : public Dump { int *vtype; // type of each vector (INT, DOUBLE) char **vformat; // format string for each vector element - int maxchoose; // size of choose arrays + int maxlocal; // size of choose and local-compute arrays int *choose; // 1 if output this atom, 0 if no double *dchoose; // value for each atom to threshhold against int nevery; int fixflag; // 1 if this dump invokes any fixes, 0 if not - int fixflag_centro; // 1 if this fix is invoked, 0 if not - int fixflag_energy; + int fixflag_epair; // 1 if this fix is invoked, 0 if not + int fixflag_centro; int fixflag_stress; - int ifix_centro; // which fix it is (0 to nfixes-1) - int ifix_energy; + int ifix_epair; // which fix it is (0 to nfixes-1) + int ifix_centro; int ifix_stress; + int computeflag; // 1 if this dump invokes any local comp, 0 if not + int computeflag_etotal; // 1 if this local comp is invoked, 0 if not + int computeflag_ke; + + // local computation arrays + + double *etotal; + double *ke; + + // private methods + int modify_param(int, char **); void write_header(int); int count(); @@ -63,6 +74,11 @@ class DumpCustom : public Dump { void header_binary(int); void header_item(int); + typedef void (DumpCustom::*FnPtrData)(int, double *); + FnPtrData write_choice; // ptr to write data functions + void write_binary(int, double *); + void write_text(int, double *); + // customize by adding a method prototype typedef void (DumpCustom::*FnPtrPack)(int); @@ -95,8 +111,10 @@ class DumpCustom : public Dump { void pack_tqx(int); void pack_tqy(int); void pack_tqz(int); + void pack_etotal(int); + void pack_ke(int); + void pack_epair(int); void pack_centro(int); - void pack_energy(int); void pack_sxx(int); void pack_syy(int); void pack_szz(int); @@ -104,10 +122,10 @@ class DumpCustom : public Dump { void pack_sxz(int); void pack_syz(int); - typedef void (DumpCustom::*FnPtrData)(int, double *); - FnPtrData write_choice; // ptr to write data functions - void write_binary(int, double *); - void write_text(int, double *); + // local computation methods + + void compute_etotal(); + void compute_ke(); }; #endif diff --git a/src/thermo.cpp b/src/thermo.cpp index 27b5890932..74fb9dd924 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -39,7 +39,7 @@ // customize by adding a keyword to this list: -// step, atoms, cpu, temp, press, pe, ke, eng, +// step, atoms, cpu, temp, press, pe, ke, etotal, // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail, erot // vol, lx, ly, lz, pxx, pyy, pzz, pxy, pxz, pyz // gke, grot (granular trans-ke and rotational-ke) @@ -48,10 +48,10 @@ // customize by adding a DEFINE to this list -#define ONE "step temp epair emol eng press" -#define MULTI "eng ke temp pe ebond eangle edihed eimp evdwl ecoul elong press" +#define ONE "step temp epair emol etotal press" +#define MULTI "etotal ke temp pe ebond eangle edihed eimp evdwl ecoul elong press" #define GRANULAR "step atoms gke grot" -#define DIPOLE "step temp epair erot eng press" +#define DIPOLE "step temp epair erot etotal press" enum {IGNORE,WARN,ERROR}; // same as write_restart.cpp enum {ONELINE,MULTILINE}; @@ -596,8 +596,8 @@ void Thermo::parse_fields(char *str) } else if (strcmp(word,"ke") == 0) { addfield("KinEng",&Thermo::compute_ke,FLOAT); tempflag = 1; - } else if (strcmp(word,"eng") == 0) { - addfield("TotEng",&Thermo::compute_eng,FLOAT); + } else if (strcmp(word,"etotal") == 0) { + addfield("TotEng",&Thermo::compute_etotal,FLOAT); tempflag = 1; } else if (strcmp(word,"evdwl") == 0) { addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT); @@ -737,10 +737,10 @@ int Thermo::compute_value(char *word, double *answer) } else if (strcmp(word,"ke") == 0) { tempvalue = temperature->compute(); compute_ke(); - } else if (strcmp(word,"eng") == 0) { + } else if (strcmp(word,"etotal") == 0) { tempvalue = temperature->compute(); fix_compute_pe(); - compute_eng(); + compute_etotal(); } else if (strcmp(word,"evdwl") == 0) compute_evdwl(); else if (strcmp(word,"ecoul") == 0) compute_ecoul(); else if (strcmp(word,"epair") == 0) compute_epair(); @@ -862,7 +862,7 @@ void Thermo::compute_ke() /* ---------------------------------------------------------------------- */ -void Thermo::compute_eng() +void Thermo::compute_etotal() { compute_pe(); double ke = 0.5 * temperature->dof * force->boltz * tempvalue; @@ -1195,7 +1195,7 @@ void Thermo::compute_pave() void Thermo::compute_eave() { - compute_eng(); + compute_etotal(); if (firststep == 0) { if (npartial_e == 0) { diff --git a/src/thermo.h b/src/thermo.h index 5de47e91b4..b2dee57df6 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -98,7 +98,7 @@ class Thermo : public LAMMPS { void compute_press(); void compute_pe(); void compute_ke(); - void compute_eng(); + void compute_etotal(); void compute_evdwl(); void compute_ecoul(); void compute_epair();