diff -Naur lammps-5Jul09/doc/Eqs/heat_flux_J.jpg lammps-6Jul09/doc/Eqs/heat_flux_J.jpg --- lammps-5Jul09/doc/Eqs/heat_flux_J.jpg 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/doc/Eqs/heat_flux_J.jpg 2009-07-02 10:38:31.000000000 -0600 @@ -0,0 +1,28 @@ +JFIFHHC  !"$"$ :<!1"AQa2q#3B$r7CRSVbu?JR)JR)JR)JR)JRmOj)/t!%Ɗ6Vs LҔ)JR)JR+JvY&U +–KKH* T@8-xjğqOv3u3[typbuUΝ'MKl'Lhcwh*IRpp O=y-F-$d!&riDD.Zv姮Pv*B:TXPI +U+ŕ(;8W6QMٛ^mj\!LOpxܞ}* ˾!C>qRyŹG-H?M[ZN¼3LA%ˍ%,pI)$}A'Bi.)ڋw7GAV ~gX.B !7$qJR/;NPoRAoJZvb -n\i-acOAI# mt$vb2n<;xQ*<``gB3whHer{H8ߥcKB]BB g)J8[YBc?ޑfK y8#pG؃G)T8)Zޞ7|3ZGv?oՆ@SaGj5i${;GY.rlx;/s'Z@+]iUkOg{7S$^6r<z7]u=j!"\"Ieevr8 r \6~v1T}?nKlsWg #̾NBR6l7z~:ɫ/_D]{Dʙ g;qf8D)y#: aAn;%[SKOTpjS3b7*:N (?Zj[m޴7Fhy)!!|7sADI k^Js:ʛ?5fhF[u=*cS[8y_wZa!'BpؔBm.%. `A|V;n04e/X Ʀ3ݒժMR_V03l`ԨM\.o2[)佄aE!ԢN>JƒԂh\-rPVut\FnZ j~2 +ZRJ7PӔpx5͏P, cًh GaDZSrCNFG$t_z$5u%,)qrܞ=W#aTF \h:Ԥi?ARкkkS;- +n0]B`zÁ1^⸆H,>@?l[Ӓ!%yBgQ>=k.)tqb&&rq c<Wnr2̆O> @[$JgNӖ۹I)op0Vyéؤ|8$%*`dx?Zul_֘SЭ͙ rBR]Z!()+R )L˜ZbF~T&[nigRSuQk+ |5U\;%YDbC4@99$8,Je##)p (`V1-Mrev +LR$ppqj +[RҢ$$$cx*,{p33r雌R6 +B#8Q8 wW+-!LdNdrG g#9vkluHSҭ̆55[uI{$) M!R +IyW:>EnpƘ{Ň9m$.}QOveθ!EzJ&GBP093Ҵ6eBC0KpϞ8Ja] z@*'>+KR +j$ R€>GlX׻Sw(mmʌ RPRqUn|qB~:.;òZS[M%IJ!,y~~nj-ea +K}"Ht惐T'=䝣EL"DvԢu%.$ 0x>٪ֽ%®L.Զ>IV9XKjZTZDuo?r=KX=%vb.v{Qܐi8SC-881މ/]mK(iFaa +q\p9'sU}bCW;$F)/%$Oa1沴 Kn lդS*YU/YX rb%ޒڴq,H)Vd]T9VF2}5ӣ[$)[Q+f:qJGG += k܅'Q{dJAs4u!) +NJ6kӶ{v#&-tGI!@ Iy>Muni)QJI&6o촔nnjs>?P]C,zAi/]$˹5(4Z.4D[ +܇-Lć-8fBFRN+{AW9zU:%rҀ 25KQeLn +$J i Ÿ\Q'yԥm*) Mon8ˋe~!&8JZHH9 JЄS q ֱZn_:L$O0ԡTǶsO֡~jw6"Zbrf5iTw +֕p`J$yW"sIjũ1U1BF@X(oj@{NL[H3&?6DnSqYu孖W>Ԅ Vj4^kCliy0 +sNH/\p K.fB4!8%kʆڱgȓ.Kn4c7m J)r + +Խ<_2}h#%[`v䒄dƛ{Jo7Y3R.l%3t$'q{7Pfr߇2诗%ˁ,6m'<: mcٖ )Y_nh{lU\nVn/|<\&"3J_ГGτT5J,V4铤ߟ5w\K#G H  +[ZZ5oi3b 6jB@I?/Ov"\e2wY@;H#;NCLۮеmBmE[ҹ!%EiA'bZ)i$VߵMj~SJr m#hc 4<;Uy`K %lJ\) q 6ϛj/ۙV훳}#oy*Fn.gZS'zy.-% Q-+'O猺QCպ:ѩqn۔R?JN/QWNiH;D"Δ<m@7 l~֮Z\Hn-_膊ARq܏Ԟ܈RCAq좒3X5DMXJJҮ (H[̭+T(JF2ϺkJS 7?PϿM'Xb=t\մV (d4ܓşC[4mx˔x,1\.)ǂTrO9F* ^hxw1sue@Khؕ%Sh%m7^g2,7gomG3TWo/O^V\!n4Ng@pm9'$sZ.\jI%d3!IBLrۉeCvxmXHݳI]G%e1B ZwZԔ9 B4eBg^iĴNvKm$'j) H# >mPnQo7;1,.41;IIZϕYR']o3>DŬFTSC%9[i @Zv4dwdȑ1 uN8.*ڐ*JR)L\']QʖRH46mv7x 4[+RB)a({4F1#cH VZRRJT`UwsqPtowS*! Z,pOyjVv BCJ#ymdޜ`q1bDiQJF}%+n|)JR)J֛o5IT1 $W".Q49!_-ӡ֠$ed{}Z j\8c)C +SLJRA궠~"Nj)̩SY);!y*_9 J|`*MLŻW#|#} @Aѵ'IJp +G2q3ksOS8B2\JNANx "Jɼ9iÄ=>"b` ˉ'9##`( +--$i{<$d}f)JR)JR)JR-@ҧǴ3 ǏSx)JRSpu}s--#ڤ6ci[\0mƈ_d8 +">m4 ز!`, zӊҶx{L1TVź#QPTq5%JR)JR xᩬׁ8ouEPukHJW)Iu#y~)JRҕHZB$>QhH6Sn2dܭ8Fi%qIZR’T0=ai}4x75&hC_Š\SK{ +@@ J!!)HRJ)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR \ No newline at end of file diff -Naur lammps-5Jul09/doc/Eqs/heat_flux_J.tex lammps-6Jul09/doc/Eqs/heat_flux_J.tex --- lammps-5Jul09/doc/Eqs/heat_flux_J.tex 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/doc/Eqs/heat_flux_J.tex 2009-07-02 10:38:31.000000000 -0600 @@ -0,0 +1,14 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ +\mathbf{J} += \sum_i e_i \mathbf{v}_i ++ \sum_{i8!1"AQ2a#Br3qUbc4R?JR)JR)JR)JR)JR)JR)JR)JRDZ$[\Hb;ON5S8+t ƸD~Tb6TY ?J)JR)JR)JW\CM-V6)J8 'T[mrgۮPDwc$q\8mlewRRI&)JR)JR)JRu u%ϚXs# !ܒaS^d wx6֘:@z< OWjM;n䈶n{6d:,˺|ԭՇrC!M )E)t2NFP4hͱ>xIhB-ClEK@RPF*HS}2)JR)JR IuDW,pYP{ۍ<$^$h ،a@P0خQD%%J 2I;}Ik%,t. G̸)%{+W;$WB&K75+6)!/HWs% z)JR)JRrh;qⵍ I I$$+Xݟ2\K !X $aXNRAc2nje$w +wPH=55% UmG<}m*F?cjPq.Gw!JN=Uz+<|i[FkvC-Ry )tZv@&!2ȋmdJ뤡}@%(U1]`yn +I{'GBa RRqa!J `+i@LB݆ˉVē10*N + ы]Dioq ;VxR1QW>qҷz.qAI_qNfj\zm_}եnٸ\Qϙ>wxdӎJ9 BIRCϼPJTRǕ'ޠ?f5"[.-GiT)Nz.t޸iG{_KC))V} d.k;2!P^qn[QRSK?uKX7oҷV%8~X +*BT'2I:ݵ/.5*SB%iS!YHHjgx+E6|[>6lhQq] ^@[PJR%DgSJ:^էf[z5Se)jϜ0dN=ԣWta{Hx;B񫂨60}&X ڬЭq;q(HHRP?!EdEy淸ڶ@|έ?2R^:S:R/C7]G! :RN<<=\r>Tk]RPoCI@He(9[} +̧ kS_)! >cb~1˷w6ۤN;F JZoIJgAY?g\j; ֛\쭶Xqm0 %Ca gbř 5Uᕩa% pssZmi׹pm6Ur̘vhgrg8x].:J˜^s閐?zVvYw%(ǛHHrA'4>UN;XVlvzٻ.QH-̅Ѵ_>TzE ]&heXq7c#85/TiMO4"Cr[enŤB*V1oZJ9=ݵ9Ʋ fB#;P󁐏Lj~n3NEQRr y`I|_8-V$y|G.Х GΤx~bU)$.#^:CoWI-{Xi-NŢkDM% Kp>R'v+Ok}]3H7  +Z[>X6n'_JRWb*\X.*" 8iK)Pbحٍ +dl%1V'WfTd@"Şښ m;G!!JBSv $>xq_=ɨ.0YX^mm۞`BqPITyHS=j[+%I)=cc$THvd';tHΩ;Je(QNsڹf= u#zRINGbFO=k5Øz;ol N0A$T,ƃS#3!.шc،QJC(l%h V jT0N⭩O8( A}ZD"`2` w<ԚRWfY`i{F-Yۊ᜷s3XbMorqz^P9;- $B6FQ\-q|Z }W CAKwP}\^ᴜ% 7vor |[P{y! eO  \9 spV6e +[je_FrlGd(I+V=Yi4ۑ-𩋕!Ʒ7$) QJJS `M9 v11 t"d% 7ԁ8ڰzjG_%FVTpy0tq83 +їˣYn-,Q.jy-KHʝrU鲨F)i<$ + +Is~ *EG +4ANs DyjEG +404_yO:]'{ JV!i +J21X qPP(m;vQ0 %!)0)JR{d1c{N%$TIDٯEUY)+).6ZҝJy)I9vCӅ:dvG"Cj3C{@ v"EM-wIw8MILJ9)NԒ88zxZ +SJi]NR DcNY丗YQ*u RԥJRTTNI$<Իm%/%z-n)kq{Br(NԤrx T)JR)JR0̨ƒamc)ZHz*FDߋjv %V +RJU봌 +$Or +4v.*e-%Sn@PhZp3h2%GөJ5q[Gt7hF^KR~d*xxQPFv䓌c.Z[N+mi)RT2pG7ZrڻuHVʷa @@BR;) +XvTcGE\RmB +SXHn>iJR)JR)JR)JR)JR)JR)JR)JR)JW \ No newline at end of file diff -Naur lammps-5Jul09/doc/Eqs/heat_flux_k.tex lammps-6Jul09/doc/Eqs/heat_flux_k.tex --- lammps-5Jul09/doc/Eqs/heat_flux_k.tex 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/doc/Eqs/heat_flux_k.tex 2009-07-02 10:38:31.000000000 -0600 @@ -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} + diff -Naur lammps-5Jul09/doc/Scripts/correlate.py lammps-6Jul09/doc/Scripts/correlate.py --- lammps-5Jul09/doc/Scripts/correlate.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/doc/Scripts/correlate.py 2009-07-02 10:38:31.000000000 -0600 @@ -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 diff -Naur lammps-5Jul09/doc/Section_commands.html lammps-6Jul09/doc/Section_commands.html --- lammps-5Jul09/doc/Section_commands.html 2009-06-22 15:17:38.000000000 -0600 +++ lammps-6Jul09/doc/Section_commands.html 2009-07-02 10:38:31.000000000 -0600 @@ -343,10 +343,10 @@

These are compute styles contributed by users, which can be used if diff -Naur lammps-5Jul09/doc/Section_commands.txt lammps-6Jul09/doc/Section_commands.txt --- lammps-5Jul09/doc/Section_commands.txt 2009-06-22 15:17:38.000000000 -0600 +++ lammps-6Jul09/doc/Section_commands.txt 2009-07-02 10:38:31.000000000 -0600 @@ -462,6 +462,7 @@ "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 -Naur lammps-5Jul09/doc/compute.html lammps-6Jul09/doc/compute.html --- lammps-5Jul09/doc/compute.html 2009-04-29 10:54:14.000000000 -0600 +++ lammps-6Jul09/doc/compute.html 2009-07-02 10:38:31.000000000 -0600 @@ -117,6 +117,7 @@

  • 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 -Naur lammps-5Jul09/doc/compute.txt lammps-6Jul09/doc/compute.txt --- lammps-5Jul09/doc/compute.txt 2009-04-29 10:54:14.000000000 -0600 +++ lammps-6Jul09/doc/compute.txt 2009-07-02 10:38:31.000000000 -0600 @@ -114,6 +114,7 @@ "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 -Naur lammps-5Jul09/doc/compute_heat_flux.html lammps-6Jul09/doc/compute_heat_flux.html --- lammps-5Jul09/doc/compute_heat_flux.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/doc/compute_heat_flux.html 2009-07-02 10:40:58.000000000 -0600 @@ -0,0 +1,152 @@ + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    compute heat/flux command +

    +

    Syntax: +

    +
    compute ID group-ID heat/flux pe-ID 
+
    +
    • ID, group-ID are documented in compute command +
    • heat/flux = style name of this compute command +
    • pe-ID = ID of a compute that calculates per-atom potential energy +
    +

    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. +

    +

    See the fix thermal/conductivity +command for details on how to compute thermal conductivity in an +alternate way, via the Muller-Plathe method. +

    +

    The compute takes a pe-ID argument which is the ID of a compute +pe/atom 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. +

    +
    +
    +
    +
    +

    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 
+
    +

    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 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: +

    +

    fix thermal/conductivity +

    +

    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 -Naur lammps-5Jul09/doc/compute_heat_flux.txt lammps-6Jul09/doc/compute_heat_flux.txt --- lammps-5Jul09/doc/compute_heat_flux.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/doc/compute_heat_flux.txt 2009-07-02 10:40:58.000000000 -0600 @@ -0,0 +1,148 @@ +"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. + +See the "fix thermal/conductivity"_fix_thermal_conductivity.html +command for details on how to compute thermal conductivity in an +alternate way, via the Muller-Plathe method. + +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:] + +"fix thermal/conductivity"_fix_thermal_conductivity.html + +[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 -Naur lammps-5Jul09/doc/fix_thermal_conductivity.html lammps-6Jul09/doc/fix_thermal_conductivity.html --- lammps-5Jul09/doc/fix_thermal_conductivity.html 2008-07-30 08:03:56.000000000 -0600 +++ lammps-6Jul09/doc/fix_thermal_conductivity.html 2009-07-02 10:40:58.000000000 -0600 @@ -52,6 +52,10 @@ Muller-Plathe method, the heat flux is imposed, and the temperature gradient is the system's response.

    +

    See the compute heat/flux command for details +on how to compute thermal conductivity in an alternate way, via the +Green-Kubo formalism. +

    The simulation box is divided into Nbin layers in the edim direction, where the layer 1 is at the low end of that dimension and the layer Nbin is at the high end. Every N steps, Nswap pairs of @@ -141,7 +145,8 @@

    Related commands:

    fix ave/spatial, fix -viscosity +viscosity, compute +heat/flux

    Default:

    diff -Naur lammps-5Jul09/doc/fix_thermal_conductivity.txt lammps-6Jul09/doc/fix_thermal_conductivity.txt --- lammps-5Jul09/doc/fix_thermal_conductivity.txt 2008-07-30 08:03:56.000000000 -0600 +++ lammps-6Jul09/doc/fix_thermal_conductivity.txt 2009-07-02 10:40:58.000000000 -0600 @@ -42,6 +42,10 @@ Muller-Plathe method, the heat flux is imposed, and the temperature gradient is the system's response. +See the "compute heat/flux"_compute_heat_flux.html command for details +on how to compute thermal conductivity in an alternate way, via the +Green-Kubo formalism. + The simulation box is divided into {Nbin} layers in the {edim} direction, where the layer 1 is at the low end of that dimension and the layer {Nbin} is at the high end. Every N steps, Nswap pairs of @@ -131,7 +135,8 @@ [Related commands:] "fix ave/spatial"_fix_ave_spatial.html, "fix -viscosity"_fix_viscosity.html +viscosity"_fix_viscosity.html, "compute +heat/flux"_compute_heat_flux.html [Default:] diff -Naur lammps-5Jul09/src/compute_heat_flux.cpp lammps-6Jul09/src/compute_heat_flux.cpp --- lammps-5Jul09/src/compute_heat_flux.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/src/compute_heat_flux.cpp 2009-07-02 10:38:31.000000000 -0600 @@ -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 -Naur lammps-5Jul09/src/compute_heat_flux.h lammps-6Jul09/src/compute_heat_flux.h --- lammps-5Jul09/src/compute_heat_flux.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jul09/src/compute_heat_flux.h 2009-07-02 10:38:31.000000000 -0600 @@ -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 -Naur lammps-5Jul09/src/style.h lammps-6Jul09/src/style.h --- lammps-5Jul09/src/style.h 2009-06-23 11:14:26.000000000 -0600 +++ lammps-6Jul09/src/style.h 2009-07-02 10:38:31.000000000 -0600 @@ -81,6 +81,7 @@ #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(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)