forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2952 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
d3cbb24f81
commit
214d346bdc
Binary file not shown.
|
After Width: | Height: | Size: 5.4 KiB |
|
|
@ -0,0 +1,14 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
|
||||
$$
|
||||
\mathbf{J}
|
||||
= \sum_i e_i \mathbf{v}_i
|
||||
+ \sum_{i<j} \left( \mathbf{f}_{ij} \cdot \mathbf{v}_j \right) \mathbf{x}_{ij}
|
||||
= \sum_i e_i \mathbf{v}_i
|
||||
+ \frac{1}{2} \sum_{i<j} \left( \mathbf{f}_{ij} \cdot \left(\mathbf{v}_i + \mathbf{v}_j \right) \right) \mathbf{x}_{ij}
|
||||
$$
|
||||
|
||||
\end{document}
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 5.0 KiB |
|
|
@ -0,0 +1,11 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
|
||||
$$
|
||||
\kappa = \frac{V}{k_B T^2} \int_0^\infty \langle J_x(0) J_x(t) \rangle \, dt
|
||||
= \frac{V}{3 k_B T^2} \int_0^\infty \langle \mathbf{J}(0) \cdot \mathbf{J}(t) \rangle \, dt
|
||||
$$
|
||||
|
||||
\end{document}
|
||||
|
||||
|
|
@ -0,0 +1,89 @@
|
|||
#!/usr/bin/python
|
||||
"""
|
||||
function:
|
||||
parse the block of thermo data in a lammps logfile and perform auto- and
|
||||
cross correlation of the specified column data. The total sum of the
|
||||
correlation is also computed which can be converted to an integral by
|
||||
multiplying by the timestep.
|
||||
output:
|
||||
standard output contains column data for the auto- & cross correlations
|
||||
plus the total sum of each. Note, only the upper triangle of the
|
||||
correlation matrix is computed.
|
||||
usage:
|
||||
correlate.py [-c col] <-c col2> <-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
|
||||
|
|
@ -343,10 +343,10 @@ description:
|
|||
</P>
|
||||
<DIV ALIGN=center><TABLE BORDER=1 >
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_centro_atom.html">centro/atom</A></TD><TD ><A HREF = "compute_cna_atom.html">cna/atom</A></TD><TD ><A HREF = "compute_coord_atom.html">coord/atom</A></TD><TD ><A HREF = "compute_damage_atom.html">damage/atom</A></TD><TD ><A HREF = "compute_displace_atom.html">displace/atom</A></TD><TD ><A HREF = "compute_erotate_asphere.html">erotate/asphere</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_erotate_sphere.html">erotate/sphere</A></TD><TD ><A HREF = "compute_group_group.html">group/group</A></TD><TD ><A HREF = "compute_ke.html">ke</A></TD><TD ><A HREF = "compute_ke_atom.html">ke/atom</A></TD><TD ><A HREF = "compute_pe.html">pe</A></TD><TD ><A HREF = "compute_pe_atom.html">pe/atom</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_pressure.html">pressure</A></TD><TD ><A HREF = "compute_reduce.html">reduce</A></TD><TD ><A HREF = "compute_reduce.html">reduce/region</A></TD><TD ><A HREF = "compute_stress_atom.html">stress/atom</A></TD><TD ><A HREF = "compute_temp.html">temp</A></TD><TD ><A HREF = "compute_temp_asphere.html">temp/asphere</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_temp_com.html">temp/com</A></TD><TD ><A HREF = "compute_temp_deform.html">temp/deform</A></TD><TD ><A HREF = "compute_temp_partial.html">temp/partial</A></TD><TD ><A HREF = "compute_temp_profile.html">temp/profile</A></TD><TD ><A HREF = "compute_temp_ramp.html">temp/ramp</A></TD><TD ><A HREF = "compute_temp_region.html">temp/region</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_temp_sphere.html">temp/sphere</A>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_erotate_sphere.html">erotate/sphere</A></TD><TD ><A HREF = "compute_group_group.html">group/group</A></TD><TD ><A HREF = "compute_heat_flux.html">heat/flux</A></TD><TD ><A HREF = "compute_ke.html">ke</A></TD><TD ><A HREF = "compute_ke_atom.html">ke/atom</A></TD><TD ><A HREF = "compute_pe.html">pe</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_pe_atom.html">pe/atom</A></TD><TD ><A HREF = "compute_pressure.html">pressure</A></TD><TD ><A HREF = "compute_reduce.html">reduce</A></TD><TD ><A HREF = "compute_reduce.html">reduce/region</A></TD><TD ><A HREF = "compute_stress_atom.html">stress/atom</A></TD><TD ><A HREF = "compute_temp.html">temp</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_temp_asphere.html">temp/asphere</A></TD><TD ><A HREF = "compute_temp_com.html">temp/com</A></TD><TD ><A HREF = "compute_temp_deform.html">temp/deform</A></TD><TD ><A HREF = "compute_temp_partial.html">temp/partial</A></TD><TD ><A HREF = "compute_temp_profile.html">temp/profile</A></TD><TD ><A HREF = "compute_temp_ramp.html">temp/ramp</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "compute_temp_region.html">temp/region</A></TD><TD ><A HREF = "compute_temp_sphere.html">temp/sphere</A>
|
||||
</TD></TR></TABLE></DIV>
|
||||
|
||||
<P>These are compute styles contributed by users, which can be used if
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -117,6 +117,7 @@ available in LAMMPS:
|
|||
<LI><A HREF = "compute_erotate_asphere.html">erotate/asphere</A> - rotational energy of aspherical particles
|
||||
<LI><A HREF = "compute_erotate_sphere.html">erotate/sphere</A> - rotational energy of spherical particles
|
||||
<LI><A HREF = "compute_group_group.html">group/group</A> - energy/force between two groups of atoms
|
||||
<LI><A HREF = "compute_heat_flux.html">heat/flux</A> - heat flux through a group of atoms
|
||||
<LI><A HREF = "compute_ke.html">ke</A> - translational kinetic energy
|
||||
<LI><A HREF = "compute_ke_atom.html">ke/atom</A> - kinetic energy for each atom
|
||||
<LI><A HREF = "compute_pe.html">pe</A> - potential energy
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -0,0 +1,147 @@
|
|||
<HTML>
|
||||
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
|
||||
</CENTER>
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
<HR>
|
||||
|
||||
<H3>compute heat/flux command
|
||||
</H3>
|
||||
<P><B>Syntax:</B>
|
||||
</P>
|
||||
<PRE>compute ID group-ID heat/flux pe-ID
|
||||
</PRE>
|
||||
<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
|
||||
<LI>heat/flux = style name of this compute command
|
||||
<LI>pe-ID = ID of a compute that calculates per-atom potential energy
|
||||
</UL>
|
||||
<P><B>Examples:</B>
|
||||
</P>
|
||||
<PRE>compute myFlux all heat/flux myPE
|
||||
</PRE>
|
||||
<P><B>Description:</B>
|
||||
</P>
|
||||
<P>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.
|
||||
</P>
|
||||
<P>The compute takes a <I>pe-ID</I> argument which is the ID of a <A HREF = "compute_pe_atom.html">compute
|
||||
pe/atom</A> 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.
|
||||
</P>
|
||||
<P>The Green-Kubo formulas relate the ensemble average of the
|
||||
auto-correlation of the heat flux J to the thermal conductivity kappa.
|
||||
</P>
|
||||
<CENTER><IMG SRC = "Eqs/heat_flux_k.jpg">
|
||||
</CENTER>
|
||||
<CENTER><IMG SRC = "Eqs/heat_flux_J.jpg">
|
||||
</CENTER>
|
||||
<P>Ei is the per-atom energy (potential and kinetic). The potential term
|
||||
is calculated by the compute <I>pe-ID</I> specified as an argument to
|
||||
the compute heat/flux command.
|
||||
</P>
|
||||
<P>IMPORTANT NOTE: The per-atom potential energy calculated by the
|
||||
<I>pe-ID</I> 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
|
||||
</P>
|
||||
<PRE>compute myPE all pe/atom pair
|
||||
</PRE>
|
||||
<P>Note that if <I>pair</I> is not listed as the last argument, it will be
|
||||
included by default, but so will other contributions such as bond,
|
||||
angle, etc.
|
||||
</P>
|
||||
<P>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.
|
||||
</P>
|
||||
<P>Here is an example of this procedure. First a LAMMPS input script for
|
||||
solid Ar is appended below. A Python script
|
||||
<A HREF = "Scripts/correlate.py">correlate.py</A> 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
|
||||
</P>
|
||||
<PRE>correlate.py flux.log -c 3 -s 200
|
||||
</PRE>
|
||||
<P>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
|
||||
<A HREF = "units.html">units</A> in LAMMPS, this is 2917703220.0 in this case. The
|
||||
final thermal conductivity value obtained is 0.25 W/mK.
|
||||
</P>
|
||||
<P><B>Output info:</B>
|
||||
</P>
|
||||
<P>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.
|
||||
</P>
|
||||
<P>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.
|
||||
</P>
|
||||
<P><B>Restrictions:</B>
|
||||
</P>
|
||||
<P>Only pairwise interactions, as defined by the pair_style command, are
|
||||
included in this calculation.
|
||||
</P>
|
||||
<P>To use this compute you must define an atom_style, such as dpd or
|
||||
granular, that communicates the velocites of ghost atoms.
|
||||
</P>
|
||||
<P><B>Related commands:</B> none
|
||||
</P>
|
||||
<P><B>Default:</B> none
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<H4>Sample LAMMPS input script
|
||||
</H4>
|
||||
<PRE>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>
|
||||
<PRE># ------------- Equilibration and thermalisation ----------------
|
||||
</PRE>
|
||||
<PRE>fix NPT all npt 70 70 10 xyz 0.0 0.0 100.0 drag 0.2
|
||||
run 8000
|
||||
unfix NPT
|
||||
</PRE>
|
||||
<PRE># --------------- Equilibration in nve -----------------
|
||||
</PRE>
|
||||
<PRE>fix NVE all nve
|
||||
run 8000
|
||||
</PRE>
|
||||
<PRE># -------------- Flux calculation in nve ---------------
|
||||
</PRE>
|
||||
<PRE>reset_timestep 0
|
||||
compute flux all heat_flux
|
||||
log flux.log
|
||||
variable J equal c_flux<B>1</B>/vol
|
||||
thermo_style custom step temp v_J
|
||||
run 100000
|
||||
</PRE>
|
||||
</HTML>
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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_{i<j} (f_ij . v_j) x_ij
|
||||
|
||||
// virial-like contribution
|
||||
// loop over neighbors of my atoms
|
||||
// require either i or j be in compute group
|
||||
|
||||
double Jv[3] = {0.0,0.0,0.0};
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
|
||||
if (j < nall) factor_coul = factor_lj = 1.0;
|
||||
else {
|
||||
factor_coul = special_coul[j/nall];
|
||||
factor_lj = special_lj[j/nall];
|
||||
j %= nall;
|
||||
}
|
||||
|
||||
if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
eng = pair->single(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);
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Reference in New Issue