diff --git a/doc/Eqs/heat_flux_J.jpg b/doc/Eqs/heat_flux_J.jpg new file mode 100644 index 0000000000..c30ae8dcf8 Binary files /dev/null and b/doc/Eqs/heat_flux_J.jpg differ diff --git a/doc/Eqs/heat_flux_J.tex b/doc/Eqs/heat_flux_J.tex new file mode 100644 index 0000000000..2c2b232058 --- /dev/null +++ b/doc/Eqs/heat_flux_J.tex @@ -0,0 +1,14 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ +\mathbf{J} += \sum_i e_i \mathbf{v}_i ++ \sum_{i <-s max_correlation_time> [logfile] +""" +import sys +import re +import array + +# parse command line + +maxCorrelationTime = 0 +cols = array.array("I") +nCols = 0 +args = sys.argv[1:] +index = 0 +while index < len(args): + arg = args[index] + index += 1 + if (arg == "-c"): + cols.append(int(args[index])-1) + nCols += 1 + index += 1 + elif (arg == "-s"): + maxCorrelationTime = int(args[index]) + index += 1 + else : + filename = arg +if (nCols < 1): raise RuntimeError, 'no data columns requested' +data = [array.array("d")] +for s in range(1,nCols) : data.append( array.array("d") ) + +# read data block from log file + +start = False +input = open(filename) +nSamples = 0 +pattern = re.compile('\d') +line = input.readline() +while line : + columns = line.split() + if (columns and pattern.match(columns[0])) : + for i in range(nCols): + data[i].append( float(columns[cols[i]]) ) + nSamples += 1 + start = True + else : + if (start) : break + line = input.readline() +print "# read :",nSamples," samples of ", nCols," data" +if( maxCorrelationTime < 1): maxCorrelationTime = int(nSamples/2); + +# correlate and integrate + +correlationPairs = [] +for i in range(0,nCols): + for j in range(i,nCols): # note only upper triangle of the correlation matrix + correlationPairs.append([i,j]) +header = "# " +for k in range(len(correlationPairs)): + i = str(correlationPairs[k][0]+1) + j = str(correlationPairs[k][1]+1) + header += " C"+i+j+" sum_C"+i+j +print header +nCorrelationPairs = len(correlationPairs) +sum = [0.0] * nCorrelationPairs +for s in range(maxCorrelationTime) : + correlation = [0.0] * nCorrelationPairs + nt = nSamples-s + for t in range(0,nt) : + for p in range(nCorrelationPairs): + i = correlationPairs[p][0] + j = correlationPairs[p][1] + correlation[p] += data[i][t]*data[j][s+t] + output = "" + for p in range(0,nCorrelationPairs): + correlation[p] /= nt + sum[p] += correlation[p] + output += str(correlation[p]) + " " + str(sum[p]) + " " + print output diff --git a/doc/Section_commands.html b/doc/Section_commands.html index b553dd271c..5a7602f8ff 100644 --- a/doc/Section_commands.html +++ b/doc/Section_commands.html @@ -343,10 +343,10 @@ description:

- - - - + + +
centro/atomcna/atomcoord/atomdamage/atomdisplace/atomerotate/asphere
erotate/spheregroup/groupkeke/atompepe/atom
pressurereducereduce/regionstress/atomtemptemp/asphere
temp/comtemp/deformtemp/partialtemp/profiletemp/ramptemp/region
temp/sphere +
erotate/spheregroup/groupheat/fluxkeke/atompe
pe/atompressurereducereduce/regionstress/atomtemp
temp/aspheretemp/comtemp/deformtemp/partialtemp/profiletemp/ramp
temp/regiontemp/sphere

These are compute styles contributed by users, which can be used if diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt index ef0ad89c7d..8ba86af7d6 100644 --- a/doc/Section_commands.txt +++ b/doc/Section_commands.txt @@ -462,6 +462,7 @@ description: "erotate/asphere"_compute_erotate_asphere.html, "erotate/sphere"_compute_erotate_sphere.html, "group/group"_compute_group_group.html, +"heat/flux"_compute_heat_flux.html, "ke"_compute_ke.html, "ke/atom"_compute_ke_atom.html, "pe"_compute_pe.html, diff --git a/doc/compute.html b/doc/compute.html index ab5641d615..1273319a88 100644 --- a/doc/compute.html +++ b/doc/compute.html @@ -117,6 +117,7 @@ available in LAMMPS:

  • erotate/asphere - rotational energy of aspherical particles
  • erotate/sphere - rotational energy of spherical particles
  • group/group - energy/force between two groups of atoms +
  • heat/flux - heat flux through a group of atoms
  • ke - translational kinetic energy
  • ke/atom - kinetic energy for each atom
  • pe - potential energy diff --git a/doc/compute.txt b/doc/compute.txt index 741eb71536..c0045acb6c 100644 --- a/doc/compute.txt +++ b/doc/compute.txt @@ -114,6 +114,7 @@ available in LAMMPS: "erotate/asphere"_compute_erotate_asphere.html - rotational energy of aspherical particles "erotate/sphere"_compute_erotate_sphere.html - rotational energy of spherical particles "group/group"_compute_group_group.html - energy/force between two groups of atoms +"heat/flux"_compute_heat_flux.html - heat flux through a group of atoms "ke"_compute_ke.html - translational kinetic energy "ke/atom"_compute_ke_atom.html - kinetic energy for each atom "pe"_compute_pe.html - potential energy diff --git a/doc/compute_heat_flux.html b/doc/compute_heat_flux.html new file mode 100644 index 0000000000..a645707a8d --- /dev/null +++ b/doc/compute_heat_flux.html @@ -0,0 +1,147 @@ + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    compute heat/flux command +

    +

    Syntax: +

    +
    compute ID group-ID heat/flux pe-ID 
+
    + +

    Examples: +

    +
    compute myFlux all heat/flux myPE 
+
    +

    Description: +

    +

    Define a computation that calculates the heat flux vector based on +interactions between atoms in the specified group. This can be used +by itself to measure the heat flux between a hot and cold reservoir of +particles or to calculate a thermal conductivity using the Green-Kubo +formalism. +

    +

    The compute takes a pe-ID argument which is the ID of a compute +pe/atom that calculates per-atom potential +energy. It should be defined for the same group used by compute +heat/flux, though LAMMPS does not check for this. +

    +

    The Green-Kubo formulas relate the ensemble average of the +auto-correlation of the heat flux J to the thermal conductivity kappa. +

    +
    +
    +
    +
    +

    Ei is the per-atom energy (potential and kinetic). The potential term +is calculated by the compute pe-ID specified as an argument to +the compute heat/flux command. +

    +

    IMPORTANT NOTE: The per-atom potential energy calculated by the +pe-ID compute should only include pairwise energy, to be consistent +with the full heat-flux calculation. Thus if any bonds, angles, etc +exist in the system, the compute should limit its calculation to only +the pair contribution. E.g. it could be defined as +

    +
    compute myPE all pe/atom pair 
+
    +

    Note that if pair is not listed as the last argument, it will be +included by default, but so will other contributions such as bond, +angle, etc. +

    +

    The heat flux J is calculated by this compute for pairwise +interactions for any I,J pair where one of the 2 atoms in is the +compute group. It can be output every so many timesteps (e.g. via the +thermo_style custom command). Then as post-processing steps, an +autocorrelation can be performed, its integral estimated, and the +Green-Kubo formula evaluated. +

    +

    Here is an example of this procedure. First a LAMMPS input script for +solid Ar is appended below. A Python script +correlate.py is also given, which calculates +the autocorrelation of the flux output in the logfile flux.log, +produced by the LAMMPS run. It is invoked as +

    +
    correlate.py flux.log -c 3 -s 200 
+
    +

    The resulting data lists the autocorrelation in column 1 and the +integral of the autocorrelation in column 2. The integral of the +correlation needs to be multiplied by V/(kB T^2) times the sample +interval and the appropriate unit conversion factors. For real +units in LAMMPS, this is 2917703220.0 in this case. The +final thermal conductivity value obtained is 0.25 W/mK. +

    +

    Output info: +

    +

    This compute calculates a vector of length 6. The 6 components are +the x, y, z components of the full heat flux, followed by the x, y, z +components of just the convective portion of the flux, which is the +energy per atom times the velocity of the atom. +

    +

    The vector values calculated by this compute are "extensive", meaning +they scale with the number of atoms in the simulation. They should be +divided by the appropriate volume to get a flux. +

    +

    Restrictions: +

    +

    Only pairwise interactions, as defined by the pair_style command, are +included in this calculation. +

    +

    To use this compute you must define an atom_style, such as dpd or +granular, that communicates the velocites of ghost atoms. +

    +

    Related commands: none +

    +

    Default: none +

    +
    + +

    Sample LAMMPS input script +

    +
    atom_style      dpd
+units 		real
+dimension	3
+boundary	p p p
+lattice 	fcc  5.376  orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
+region  	box block 0 4 0 4 0 4
+create_box 	1 box
+create_atoms 	1 box
+mass 		1 39.948
+pair_style	lj/cut 13.0
+pair_coeff	* * 0.2381 3.405
+group 		every region box
+velocity 	all create 70 102486 mom yes rot yes dist gaussian
+timestep 	4.0
+thermo	        10 
+
    +
    # ------------- Equilibration and thermalisation ---------------- 
+
    +
    fix 		NPT all npt 70 70 10 xyz 0.0 0.0 100.0 drag 0.2
+run 		8000
+unfix           NPT 
+
    +
    # --------------- Equilibration in nve ----------------- 
+
    +
    fix 		NVE all nve
+run 		8000 
+
    +
    # -------------- Flux calculation in nve --------------- 
+
    +
    reset_timestep 0
+compute 	flux all heat_flux
+log     	flux.log
+variable        J  equal c_flux1/vol
+thermo_style 	custom step temp v_J 
+run 	        100000 
+
    + diff --git a/doc/compute_heat_flux.txt b/doc/compute_heat_flux.txt new file mode 100644 index 0000000000..1cd1727ff1 --- /dev/null +++ b/doc/compute_heat_flux.txt @@ -0,0 +1,142 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute heat/flux command :h3 + +[Syntax:] + +compute ID group-ID heat/flux pe-ID :pre + +ID, group-ID are documented in "compute"_compute.html command +heat/flux = style name of this compute command +pe-ID = ID of a compute that calculates per-atom potential energy :ul + +[Examples:] + +compute myFlux all heat/flux myPE :pre + +[Description:] + +Define a computation that calculates the heat flux vector based on +interactions between atoms in the specified group. This can be used +by itself to measure the heat flux between a hot and cold reservoir of +particles or to calculate a thermal conductivity using the Green-Kubo +formalism. + +The compute takes a {pe-ID} argument which is the ID of a "compute +pe/atom"_compute_pe_atom.html that calculates per-atom potential +energy. Normally, it should be defined for the same group used by +compute heat/flux, though LAMMPS does not check for this. + +The Green-Kubo formulas relate the ensemble average of the +auto-correlation of the heat flux J to the thermal conductivity kappa. + +:c,image(Eqs/heat_flux_k.jpg) + +:c,image(Eqs/heat_flux_J.jpg) + +Ei is the per-atom energy (potential and kinetic). The potential +portion is calculated by the compute {pe-ID} specified as an argument +to the compute heat/flux command. + +IMPORTANT NOTE: The per-atom potential energy calculated by the +{pe-ID} compute should only include pairwise energy, to be consistent +with the second virial-like term in the formula for J. Thus if any +bonds, angles, etc exist in the system, the compute should limit its +calculation to only the pair contribution. E.g. it could be defined +as follows. Note that if {pair} is not listed as the last argument, +it will be included by default, but so will other contributions such +as bond, angle, etc. + +compute myPE all pe/atom pair :pre + +The second term of the heat flux equation for J is calculated by +compute heat/flux for pairwise interactions for any I,J pair where one +of the 2 atoms in is the compute group. It can be output every so +many timesteps (e.g. via the thermo_style custom command). Then as +post-processing steps, an autocorrelation can be performed, its +integral estimated, and the Green-Kubo formula evaluated. + +Here is an example of this procedure. First a LAMMPS input script for +solid Ar is appended below. A Python script +"correlate.py"_Scripts/correlate.py is also given, which calculates +the autocorrelation of the flux output in the logfile flux.log, +produced by the LAMMPS run. It is invoked as + +correlate.py flux.log -c 3 -s 200 :pre + +The resulting data lists the autocorrelation in column 1 and the +integral of the autocorrelation in column 2. The integral of the +correlation needs to be multiplied by V/(kB T^2) times the sample +interval and the appropriate unit conversion factors. For real +"units"_units.html in LAMMPS, this is 2917703220.0 in this case. The +final thermal conductivity value obtained is 0.25 W/mK. + +[Output info:] + +This compute calculates a vector of length 6. The 6 components are +the x, y, z components of the full heat flux, followed by the x, y, z +components of just the convective portion of the flux, which is the +energy per atom times the velocity of the atom. + +The vector values calculated by this compute are "extensive", meaning +they scale with the number of atoms in the simulation. They should be +divided by the appropriate volume to get a flux. + +[Restrictions:] + +Only pairwise interactions, as defined by the pair_style command, are +included in this calculation. + +To use this compute you must define an atom_style, such as dpd or +granular, that communicates the velocites of ghost atoms. + +[Related commands:] none + +[Default:] none + +:line + +Sample LAMMPS input script :h4 + +atom_style dpd +units real +dimension 3 +boundary p p p +lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 +region box block 0 4 0 4 0 4 +create_box 1 box +create_atoms 1 box +mass 1 39.948 +pair_style lj/cut 13.0 +pair_coeff * * 0.2381 3.405 +group every region box +velocity all create 70 102486 mom yes rot yes dist gaussian +timestep 4.0 +thermo 10 :pre + +# ------------- Equilibration and thermalisation ---------------- :pre + +fix NPT all npt 70 70 10 xyz 0.0 0.0 100.0 drag 0.2 +run 8000 +unfix NPT :pre + +# --------------- Equilibration in nve ----------------- :pre + +fix NVE all nve +run 8000 :pre + +# -------------- Flux calculation in nve --------------- :pre + +reset_timestep 0 +compute flux all heat_flux +log flux.log +variable J equal c_flux[1]/vol +thermo_style custom step temp v_J +run 100000 :pre + diff --git a/src/compute_heat_flux.cpp b/src/compute_heat_flux.cpp new file mode 100644 index 0000000000..320cab75a2 --- /dev/null +++ b/src/compute_heat_flux.cpp @@ -0,0 +1,225 @@ +/* ---------------------------------------------------------------------- + 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-93AL85000 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: Reese Jones (Sandia) + Philip Howell (Siemens) + Vikas Varsney (Air Force Research Laboratory) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "compute_heat_flux.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "force.h" +#include "pair.h" +#include "modify.h" +#include "group.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM}; + +/* ---------------------------------------------------------------------- */ + +ComputeHeatFlux::ComputeHeatFlux(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all("Illegal compute heat/flux command"); + + vector_flag = 1; + size_vector = 6; + extvector = 1; + + // store pe/atom ID used by heat flux computation + // insure it is valid for pe/atom computation + + int n = strlen(arg[3]) + 1; + id_atomPE = new char[n]; + strcpy(id_atomPE,arg[3]); + + int icompute = modify->find_compute(id_atomPE); + if (icompute < 0) error->all("Could not find compute heat/flux pe/atom ID"); + if (modify->compute[icompute]->peatomflag == 0) + error->all("Compute heat/flux compute ID does not compute pe/atom"); + + vector = new double[6]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeHeatFlux::~ComputeHeatFlux() +{ + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHeatFlux::init() +{ + // error checks + + if (atom->avec->ghost_velocity == 0) + error->all("Compute heat/flux requires ghost atoms store velocity"); + + if (force->pair == NULL || force->pair->single_enable == 0) + error->all("Pair style does not support compute heat/flux"); + + int icompute = modify->find_compute(id_atomPE); + if (icompute < 0) + error->all("Pe/atom ID for compute heat/flux does not exist"); + atomPE = modify->compute[icompute]; + + pair = force->pair; + cutsq = force->pair->cutsq; + + // need an occasional half neighbor list + + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHeatFlux::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHeatFlux::compute_vector() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double rsq,eng,fpair,factor_coul,factor_lj,factor; + double fdotv,massone,ke,pe; + int *ilist,*jlist,*numneigh,**firstneigh; + + invoked_vector = update->ntimestep; + + double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + + neighbor->build_one(list->index); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // heat flux J = \sum_i e_i v_i + \sum_{isingle(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); + + if (newton_pair || j < nlocal) factor = 1.0; + else factor = 0.5; + + // symmetrize velocities + + double vx = 0.5*(v[i][0]+v[j][0]); + double vy = 0.5*(v[i][1]+v[j][1]); + double vz = 0.5*(v[i][2]+v[j][2]); + fdotv = factor * fpair * (delx*vx + dely*vy + delz*vz); + + Jv[0] += fdotv*delx; + Jv[1] += fdotv*dely; + Jv[2] += fdotv*delz; + } + } + } + + // energy convection contribution + // uses per-atom potential energy + + if (!(atomPE->invoked_flag & INVOKED_PERATOM)) { + atomPE->compute_peratom(); + atomPE->invoked_flag |= INVOKED_PERATOM; + } + + double *mass = atom->mass; + double *rmass = atom->rmass; + double mvv2e = force->mvv2e; + + double Jc[3] = {0.0,0.0,0.0}; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + massone = (rmass) ? rmass[i] : mass[type[i]]; + ke = mvv2e * 0.5 * massone * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + pe = atomPE->scalar_atom[i]; + eng = pe + ke; + Jc[0] += v[i][0]*eng; + Jc[1] += v[i][1]*eng; + Jc[2] += v[i][2]*eng; + } + } + + // total flux + + double data[6] = {Jv[0]+Jc[0],Jv[1]+Jc[1],Jv[2]+Jc[2], + Jc[0],Jc[1],Jc[2]}; + MPI_Allreduce(data,vector,6,MPI_DOUBLE,MPI_SUM,world); +} diff --git a/src/compute_heat_flux.h b/src/compute_heat_flux.h new file mode 100644 index 0000000000..54aac14d79 --- /dev/null +++ b/src/compute_heat_flux.h @@ -0,0 +1,39 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_HEATFLUX_H +#define COMPUTE_HEATFLUX_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeHeatFlux : public Compute { + public: + ComputeHeatFlux(class LAMMPS *, int, char **); + ~ComputeHeatFlux(); + void init(); + void init_list(int, class NeighList *); + void compute_vector(); + + private: + double **cutsq; + class Pair *pair; + class NeighList *list; + class Compute *atomPE; + char *id_atomPE; +}; + +} + +#endif diff --git a/src/style.h b/src/style.h index c6c6560411..1a5afd7237 100644 --- a/src/style.h +++ b/src/style.h @@ -81,6 +81,7 @@ CommandStyle(write_restart,WriteRestart) #include "compute_coord_atom.h" #include "compute_displace_atom.h" #include "compute_group_group.h" +#include "compute_heat_flux.h" #include "compute_ke.h" #include "compute_ke_atom.h" #include "compute_pe.h" @@ -106,6 +107,7 @@ ComputeStyle(cna/atom,ComputeCNAAtom) ComputeStyle(coord/atom,ComputeCoordAtom) ComputeStyle(displace/atom,ComputeDisplaceAtom) ComputeStyle(group/group,ComputeGroupGroup) +ComputeStyle(heat/flux,ComputeHeatFlux) ComputeStyle(ke,ComputeKE) ComputeStyle(ke/atom,ComputeKEAtom) ComputeStyle(pe,ComputePE)