Notes
284c1f2d04aclammps-icms

Notes
View Options

			Documentation and other notes on LAMMPS

			--------------------------------------------
			valgrind run:

			valgrind --leak-check=yes --show-reachable=yes --num-callers=10 lmp_g++ < in.*

			--------------------------------------------
			10 Jan 2007
			size of LAMMPS

			85K lines, 63K are in style which is 75%
			rest is "core" of LAMMPS

			--------------------------------------------
			1 Dec 2006
			features wish list

			Pieter's new stuff:
			-fix: Cross-linking scheme (including updating angles, torsions, etc.)
			-pair_style: Ewald summation for various types of LJ and Buckingham
			-pair_style: interaction between spheres made of solid LJ-ium
			-box: triclinic PBCs
			-fix: SLLOD algorithm for NEMD (requries triclinic PBCs)
			-fix: 1D density and velocity profiles

			speed optimizations: tune individual kernels (pair potentials,
			neighbor list construction), learn from GROMACS how they are so fast
			(assembler, better lo-level data structures), implement kernels for
			multi-core nodes or vector (SIMD-style) nodes

			parallel scalability enhancements: implement some of IBM or DE Shaw
			algorithmic ideas to enable better scaling for small problems
			on large numbers of processors

			triclinic boundary conditions: non-orthogonal simulation boxes,
			ability to model crystalline lattices with those asymmetries, do bulk
			shear viscosity modeling

			point-dipole force fields: short-range is done, but
			long-range (Ewald and PPPM) is challenging

			aspherical particle models: ellipsoidal, multiple glued spheres, DEM -
			working on with Schunk/Grest groups (CRADA, NNEDC, LDRD)

			colliodal systems mixed with small-particle solvents: new potentials,
			atom types, long-range interactions, efficiency issues for mixed-size
			systems - working on with Schunk/Grest groups

			atom-fluid hybrid models: how to model particles in a background fluid
			(i.e. nanoparticle suspensions), have to couple to an efficient fluid
			solver - working on with Schunk/Grest groups

			coupling LAMMPS to GAMESS (quantum) or FE (continuum) codes: work with
			coarse-grain continuum particles (PeriDynamics) and directly coupling
			to FE meshes (various methods exist, Lehoucq MICS project, Mark
			Robbins (JHU) work)

			other coarse-grain particle models that could be supported by LAMMPS:
			DSMC, SPH

			add modified EAM (MEAM) potential: working with Greg Wagner

			add universal REAXX potential: Aidan working with Bill Goddard group

			add charge-equilibration ability: GRASP feature that Aidan wants
			to port to LAMMPS

			polarizable force field: several variants exist

			AI-REBO force field for C nanotubes and their interaction with CH
			surroundings: worked on by Ase Henry (summer student) but not finished


			--------------------------------------------
			10 Oct 2006
			new WWW site and source-code repository

			For changing just doc and/or www pages:

			a) check all changes in doc,www into the repository
			b) run lmpwww.py

			For fixing a bug w/ a new tarball:

			a) check all bug fix src,doc,tools,etc changes into the repository
			make SURE that "make package-check" shows all package dirs in sync
			repository contains NO installed packages and NO package files in src dir,
			so changes should ALWAYS only go into package dirs
			NEVER do "svn add" for a package file that is in the top-level src dir
			b) edit www/bug.txt to describe the new patch and run txt2html on it
			c) run lmptar.py to create the new tarball which will push it onto the
			WWW server, and put patch files in my dir
			d) svn add the 2 patch files to the repository
			e) check in www/bug.txt/html and the 2 patch files from www/patch
			into the repository
			f) run lmpwww.py to update the WWW with the bug fix

			useful SVN commands (from any level within lammps dir)

			svn help command get help on a command and its options
			svn -qu status * files are different in repository
			M files are different on my machine
			-u show update info, -q = print as little as possible
			svn add file.cpp add a file, so repository will know about it
			svn delete file.cpp delete a file, so repository will know it's gone
			svn diff file.cpp local diff of my changed file versus version I checked out
			svn -r N diff file.cpp diff of my file versus current repos (get N from stat)
			svn update grab newer files from repository and overwrite mine
			svn ci -m "" check in my changed files into repository with NULL comment

			WWW interface to LAMMP repository:
			https://development.sandia.gov/viewvc/lammps

			--------------------------------------------
			21 Sept 2006
			how setflag works in pair potentials

			setflag only used for i >= j (half of all i,j)

			what top-level pair.cpp does with setflag:
			when init(), checks that setflag i,i is set for all atom types
			leaves it up to init_one(i,j) to check setflag i,j values

			most pair styles:
			setflag is created in allocate() and all i,j set to 0
			happens when 1st coeff command is called by data file or input script
			when pair_style command is re-invoked (allocation already done):
			if setflag i,j is set, cut_i,j is reset to new global cut
			when coeff() is called with range of i,j:
			setflag i,j is set for each i,j in range
			when init_one is called for i,j:
			if setflag i,j = 0, mixing is performed
			when writing restart file:
			loop over all i,j, only write coeffs if setflag i,j is set

			some pair styles do not allow mixing (have to specify i,j explicitly):
			buck, buck/coul/cut, buck/coul/long, dpd, morse, peri, table

			table pair style:
			destroys setflag when pair_style is re-invoked:
			this forces user to reset all coeffs

			granular pair styles:
			coeff() is illegal to call, so when pair_style is called
			it sets all i,j setflag values for all atom types
			in hybrid mode, this is not just types gran is used for
			since it has no idea of these since coeff() is not called
			but it should be OK to set setflag for all i,j

			pair eam:
			only sets i,i values of setflag in coeff
			also sets mass for i,i calls to coeff
			it will mix i,j, but isn't done until init_style,
			after all i,j calls to init_one(i,j)
			so init_one for i,j call tests if setflag i,i and j,j are set
			else an error since won't be able to mix later

			potentials that requires one call to coeff with * * and allows
			some elements to be masked out:
			eam/alloy, eam/fs, sw, tersoff, airebo
			coeff() also requires a mapping list of atom types to elements
			clears setflag every time coeff() is called
			sets setflag only if I,J are both non-NULL in mapping list
			also sets mass if I = J

			pair hybrid:
			destroys setflag when pair_style is re-invoked:
			this forces user to reset all coeffs
			if coeff() called with one_coeff sub-style (requires * * and uses
			type maps), then unsets all I,J previously set to that sub-style
			this is so a 2nd * * eam/alloy call will wipe out settings from 1st one
			coeff() calls sub-style with same range of I.J args it was called with
			coeff() loops over I,J it called with and only sets setflag I,J
			if the sub-style set I,J or if sub-styly = "none" (NULL)
			sets setflag for NULL, so that "pair_coeff I J none" can be used
			and all setflag values will be set by entirety of pair_coeff commands

			--------------------------------------------
			11 Sept 2006
			reworking of neighbor lists in LAMMPS to allow 1/2 and full together

			neighbor.cpp and its children (half, full, etc) all got changed
			half routines always do build in half data structs
			full routines always do build in full data structs
			neigh_half, neigh_full flag whether half/full are ever built
			knows which (or both) to build every N steps
			has setup/cleanup to allow half/full to be called even if not init()
			e.g. by delete_atoms
			pair routines now flag what kind of neigh list they want
			neigh_half_every, neigh_full_every
			never needs to call neigh build since neigh will always build them
			pair hybrid might set both
			it also accesses correct neigh list when building sub-list
			pair peri sets neither
			pair airebo sets full
			pair sw sets full
			fix routines can flag what kind of neigh list they want
			neigh_half_every, neigh_half_once, neigh_full_once, neigh_full_every
			then calls neigh build if needed
			fix centro wants full once
			fix energy and fix stress and fix rdf want half once
			fix orient/fcc wants full every
			delete_atoms needs a half list
			so calls build_setup and cleanup (if neighbor doesn't do half)
			also calls build_half directly

			--------------------------------------------
			1 Sept 2006
			discussion with Aidan about LAMMPS + GRASP

			Hi Steve,

			Good analysis. I put some comments below. I like the idea in principle,
			because there are lots of productivity gains. It was good for me to write my
			own MD code, because I learned a lot, I had a lot of flexibility to try new
			things, and it got me a lot of visibility. All of these things are still
			arguments for continuing GRASP, but the benefits are diminishing.

			The thing that worries me is the transition process. The worst case
			scenario is that I might end up supporting two codes.

			Also, I don't want to upset existing GRASP users.

			I need to figure out a transition path that will mimimize the risk of
			getting stuck in a big hole, and it can't get in the way of my existing
			near-term commitments. The two big items are:

			- Implement a pressure calculation for ReaxFF
			- Shock simulation capability:
			* Threebody force field (in GRASP, not in LAMMPS)
			* Special initial/boundary conditions (not in GRASP, probably in LAMMPS)
			* On-the-fly analysis (not in either)

			I am thinking that for now, I will continue to use GRASP for ReaxFF and I
			will continue to primarily suppport GRASP. But I would like to get the
			threebody force field into LAMMPS. That should be relatively easy, and it
			will be a good test case. It will also take care of those people at Sandia
			who currently use GRASP for Stillinger-Weber Silicon.

			Depending on how that goes, we can then consider porting Tersoff, Charge
			Equilibration and ReaxFF, at which point GRASP can be moth-balled.

			If we port ReaxFF to LAMMPS, there might be a lot of things we want to do
			differently (including ditching the FORTRAN code and making ReaxFF open
			source). So it is good to know all the issues first.

			So for now, let's not publicize a possible merger/buy-out, but just say that
			we will port the threebody force field from GRASP to LAMMPS.

			Aidan

			> Aidan - I thought some more about what it would mean to have all of
			> LAMMPS + GRASP in one code. Did you think of more features in
			> GRASP that would be important to have in the combined code?
			>
			> My thoughts are below.
			>
			> Steve
			>
			> ----------------------
			>
			> LAMMPS + GRASP
			>
			> 1) Main benefits
			>
			> One code to maintain, advertise, distribute, attract users to.
			> One code with all features so we and others can do more problems.
			> One code with all these features becomes more unique, valuable.

			Agreed.

			>
			> 2) Political issues
			>
			> What to name the combined code?
			> LAMMPS, LAMMPS/GRASP, some new name, etc
			> My vote would be just LAMMPS since it is now a reasonably well-known
			> MD code, but that is a biased view.
			> It could be bad for you and GRASP if it loses its separate identity.
			> We need to somehow make it clear the new code is
			> significantly different/enhanced beyond the original LAMMPS.

			So here are the features that would make LAMMPS "new and improved"

			Non-optional stuff (things I must support):
			-ReaxFF
			-Charge Equilibration (required for ReaxFF, but indepenendent of it)
			-Threebody: requires complete neighbor list
			-Tersoff (each ij pair interaction is modified by a function that is a sum
			of ijk interactions, j != k)

			Optional stuff:
			-Triclinic pbc's (or is it in LAMMPS now?)
			I think a lot of the special features in GRASP (e.g. initializing atom
			positions based on a density profile, selective thermostatting) are in
			LAMMPS or could be added fairly easily using existing framework.

			>
			> Who "owns" the new code?
			> Simple answer: all of the developers in 1400 (you, Paul, me).
			> Make the new code's WWW page not be Plimpton-centric.
			> I.e. on the home page and "Authors"
			> page we list all of our names (me, you, Paul) as having created
			> the code and as contacts with email addresses, etc.
			> We all could monitor the mail list, etc.

			This is still tricky, for two usual reasons: credit and work. First,
			everyone likes to get recognized for individual accomplishments, and this is
			harder to do if the code has a collective ownership. Second, just
			maintaining the code will require on-going work, and there will have to be
			some way to divide this up. I am happy to field questions and fix bugs
			related to the GRASP force fields, but then what happens when someone wants
			to use charge equilibration with CHARMM?

			> This would also be the story we tell to management. I.e. We have combined
			> the codes to make one new MD code with all the features of both.
			> It is now no longer Plimpton's code, but 1400s.
			>
			> 3) Technical issues:
			>
			> Hooking in new pair potentials (SW, Tersoff, REAXX) to LAMMPS is
			> potentially easy, if they can be massaged into the following interface
			> by which LAMMPS uses all its pair potentials:
			>
			> input commands:
			> pair_style tersoff arg1 arg2 ...
			> args are any global settings, e.g. cutoff
			> pair_coeff i j file 1 2 3 4 (similar to EAM potentials)
			> i,j are atom types, can be , to specify all types
			> file could be one or more files that specify the force field
			> 1,2,3,4 are a "mapping" that maps LAMMPS atom types to elements
			> in the force field, e.g. Si for Tersoff, or something else in REAXX
			>
			ReadInteractions()

			> setup routines:
			> one-time setup of potential
			> i.e. read files, do pre-computation
			>
			SetupInteractions()

			> init()
			> setup the potential each time before a run
			> usually just checks that all pair_coef settings have been made
			> and figures out cutoffs for each atom type pair
			>
			SetupNeighbor()

			> compute()
			> called every timestep
			> has a neighbor list (half, full, whatever) that it loops over
			> computes force on each atom, energy (optional), virial (optional)
			ApplyForce()

			> single()
			> called by diagnostic functions that compute per-atom energy, stress
			> compute the I-J interaction (force, energy) for 2 atoms
			> probably ignore this for complicated many-body potentials
			Good idea, currently per-atom energy is done (or not) in the energy
			calculation.

			> restart()
			> called by restart files
			> how to read/write the potential parameters to a restart file
			> probably ignore this for a complicated potential
			GRASP does not have this

			>
			> 4) Distributing REAXX:
			>
			> My vote would be to hook LAMMPS to REAXX via library calls and to
			> treat it in LAMMPS as an optional "package". This is what we do
			> already for the external package POEMS (rigid body dynamics). Thus
			> REAXX could be distributed with LAMMPS in a separate directory (or
			> not, if the CalTech people don't want it to be).
			>
			> So if someone wants to use REAXX, they first build REAXX (with its own
			> Makefile) as a library. Then they include the LAMMPS REAXX package
			> and build LAMMPS, linking to the library. There would be a thin
			> pair_reaxx.cpp class in LAMMPS that made calls to the REAXX library.
			>
			> I think this would also avoid any issues with REAXX being Fortran,
			> since once it is built as a library, LAMMPS doesn't care what language
			> it is written in - it just has a callable interface.

			Sounds good.

			--------------------------------------------
			21 Jul 2006
			special bond treatment in LAMMPS vs other packages

			see email exchange with Paul and Naveen and Marcus in NOTES/special_bonds

			answer: CHAMM always applies 1-2,1-3,1-4 weighting based on bond topology,
			not on which angle, dihedrals happen to be defined

			for 1-4, it actually applies a weighting that depends on the types of
			the 2 atoms (and possibly on what comes in between?, though this
			doesn't make sense)

			it has no ambiguity in 5- or 6-atom rings, since the 2 atoms are either
			a 1-3 or 1-4 and the correct weighting gets applied directly

			in LAMMPS, we apply the pairwise weighting based on bond topology as
			well, however we apply the speical weighting as part of the dihedral
			interaction - we read in the special sigma/epsilon as part of pair
			coeffs (and it gets mixed there), but it gets computed in the dihedral
			with an additional weighting term that is an added dihedral coeff

			the CHARMM pre-proc takes care to do this correctly in special cases:
			a) defines a K = 0 dihedral if there needs to be a 1-4 LJ term
			but there is no dihedral defined
			b) 5- and 6-atom rings - e.g. in the 6-atom case, the 1-4 atoms appear
			in 2 dihedrals, but the LJ term should only apppear once, so the
			dihedral weight is set to 1/2 its value

			--------------------------------------------
			20 Jul 2006
			neighbor binning details

			we use global bin numbering scheme
			this lets us fit more, smaller bins than if tried to have integer
			# of bins in each proc's sub-domain

			global bin numbering scheme in each dim:
			0 = 0.0 to binsize
			1 = binsize to 2*binsize
			nbin-1 = prd-binsize to binsize
			nbin = prd to prd+binsize, etc
			-1 = -binsize to 0.0, etc

			mbinlo in each dim for each proc must be leftmost bin that a owned
			or ghost particle will ever be in at time neigh lists are formed
			doesn't matter that mbinlo or bin #'s go negative since
			each proc will treat mbinlo as a local zero
			and difference its bins from mbinlo to yield a positive integer index

			then I add 1 to extend mbinlo and mbinhi by 1 to cover stencil bins
			this is to insure a stencil bin index never gets generated that would
			be outside mbinlo/hi when applied to the bin of an owned atom
			not sure why delta = 1 is the correct value
			if move to center-of-interaction algorithm for communicating 1/2
			the cutoff distance, then will generated stencil bins around
			ghost atom bins and will likely need to extend mbinlo/hi by more
			than 1
			best way would be to check that stencil is fully contained in mbinxlo/hi
			after compute it, for any owned (or ghost) atom

			choose bin size to be roughly 1/2 of cutneigh, b/c it appears to be
			roughly optimal in terms of stencil size
			if bins are much smaller, then bins have few or no atoms and loop
			over stencil for each atom is too costly
			if bins are much bigger, then bins have lots of atoms and are testing
			on many atoms that are too far away when form neighbor list
			note that code should work no matter what bin size is chosen since
			next(xyz) and stencil extend to accommodate any bin size
			this is why even if box size is tiny and cutoff is huge the code
			still forms correct neighbor lists
			in this case, the bin size has to be >= 1, and binsize will be < cut/2
			but next(xyz) will then be large and the stencil will extend for many
			bins to find all the neighbors

			MUST be exactly integer # of bins in each box dimension
			this is b/c of the way bins are used in the neighbor-list formation
			i.e. loop over 13 or 26 bins, check atoms to the right and above, etc
			bin usage requires that 2 neighboring procs always think the bin
			boundary is in same place
			this is obviously the case for interior bins since a global bin system
			is used
			also needs to be the case for procs at end of box that are neighbors
			thru PBC
			if weren't the case, here is example that breaks code for newton bin
			if an integer # of bins is not enforced:
			imagine 2 nearby atoms, straddling PBC in x, with the righmost one
			having a smaller z value
			the proc on the left end will not consider them neighbors, b/c
			the ghost atom is in a bin to the left, and thus not looped over
			the proc on the right end will not consider them neighbors if the 2
			atoms are in the same bin, b/c the rightmost atom has a smaller
			z coord
			the proc on the right end will consider them neighbors IF the
			ghost atom (rightmost) is in a different bin
			this implies that the bin boundary must be at the right edge
			of the box, to be consistent with the bin boundary at the left edge
			of the box

			this is also why there is a careful test in coord2bin to insure
			atoms get put in the correct bins even when there might be round-off errors
			this test is based on there being exactly an integer # of bins that
			fit in the periodic box
			any atom with coord >= boxhi must be a ghost atom
			any atom with coord < boxlo must be a ghost atom
			so put those kinds of atoms directly in ghost bins
			put atoms in between (owned atoms) directly in non-ghost bins

			--------------------------------------------
			12 Jul 2006
			PPPM variable details

			set_grid() sets nx_pppm, ny_pppm, nz_pppm = size of PPPM grid
			includes slab size, since zprd is scaled by slab size

			nxlo_in,nxhi_in
			nylo_in,nyhi_in
			nzlo_in,zyhi_in = extent of PPPM grid that I own
			z values done for smaller zprd grid, not zprd_slab

			nxlo_out,nxhi_out
			nylo_out,nyhi_out
			nzlo_out,zyhi_out = extent of PPPM grid that my particles contribute charge to
			z values done for smaller zprd grid, not zprd_slab

			upper end +z procs have their nzhi_in and nzhi_out reset to nz_pppm-1

			nxlo_ghost,nxhi_ghost
			nylo_ghost,nyhi_ghost
			nzlo_ghost,nzhi_ghost = # of planes I recv from each dir
			to overlay domain I own

			-z proc sends to +z proc, but not vice versa

			--------------------------------------------
			5 Jun 2006
			output rules for LAMMPS

			3 kinds of output:
			thermo = always 1st/last/every-N step
			never repeat an output twice in same run
			dump = always 1st/last/every-N but never repeat a step in file
			restart = not-1st/last/every-N
			can also do anytime via write_restart command in input script
			special cases to test these outputs for
			0-step run, non-N run, N-run, N-exit min, early-exit min, 0-iter min
			thermo N with N = 0 or N > 0
			using reset_timestep between runs/min to change timesteps
			never double print thermo or dump or restart

			--------------------------------------------
			17 May 2006
			humorous email requests

			Also, When I tried downloading LAMMPS simulator software it did in
			.TAR format. I'm not able to find software to untar the tar files. Can
			you please suggest any software for that?

			--------------------------------------------
			14 May 2006
			temperature usage

			places in LAMMPS that a Temperature class is defined and how it is used:

			force templist
			defines "default" when code starts up
			any input script "temperature" command creates new one, added to templist
			any fix command that calls add_temp uses flag = 1 to allow for overwrite
			of existing one
			this is needed if fix in unfix, then recreated, since temp ID will be same
			ditto if fix is redefined
			none of entries in templist are ever deleted
			don't want to do this since don't know all classes which are storing
			ptrs to that temp (e.g. thermo, other fixes, etc)
			force does init of all in templist every run, deletes them at end

			fix nph, npt, nvt, temp/rescale
			each define a new temperature in force::templist when created via add_temp
			nph, npt = temperature fix-ID all full
			nvt = temperature fix-ID fix-group-ID full
			temp/rescale = temperature fix-ID fix-group-ID full/region
			full/region depends on if fix parameter region was used
			if user later specifies fix_modify temp ID then internal temperature ptr
			is changed to point to new one from templist
			init/delete on temp never called, let force do it

			velocity
			only used for duration of single command
			creates/deletes its own (temporarily) as temperature temp vel-group full
			or uses one from force->templist if user specifies
			calls init() on whichever one it is
			this means fixes that subtract dof need to already be defined
			e.g. fix shake or fix rigid

			thermo
			uses force->templist[0] = default by default
			user can override this with thermo_modify temp ID
			creates/deletes/init on temperature never called, let force do it
			only calls t->compute() and t->dof()
			used for fields: temp, ke, press, total E
			is user error if assigns T to thermo that is inconsistent with global P
			can also do thermo print out of any T by using keyword t_ID
			this just calls t->compute for that temp in force->templist

			pressure
			just gets passed in by whoever calls press
			pressure does not define it's own temperature
			pressure is defined/used by:
			force - the master definition
			fix npt, nph (pass their own or user-defined T to pressure)
			thermo (passes default T to pressure)

			--------------------------------------------
			14 May 2006
			dump features

			y = supported
			- = ignored (not implemented, but not an error)
			err = throws an error

			style atom bond custom dcd xtc xyz

			format y y y - - y
			scale y - - - - -
			image y - - - - -
			flush y y y y - y
			every y y y err y y
			region - - y - - -
			thresh - - y - - -

			.bin y - y - - -
			.gz y y y - - y
			per-proc y y y - - -
			per-snap y y y - - y

			only only
			group any any any all all any
			sorted no no no y y y (if group = all)

			--------------------------------------------
			2 May 2006
			how to use gprof

			compile/link with -g -pg -O (e.g. in Makefile.g++)
			run a sample input script (in.lj) which produces gmon.out.* files
			gprof lmp_g++ gmon.out.* > gprof.report

			--------------------------------------------
			17 Apr 2006
			lattice usage

			xlattice, ylattice, zlattice only used as distane conversion
			factors from box <-> lattice units by a few commands:

			domain.cpp - creates settings and maps lattice <-> box with them
			displace_atoms.cpp
			fix_indent.cpp
			region.cpp
			temperature.cpp
			velocity.cpp
			create_atoms.cpp - uses implicitly since it calls box2lattice and lattice2box

			--------------------------------------------
			17 Apr 2006
			discussion w/ Gary Grest about viscosity drag term

			Steve, I am not sure how most people think of the drag term. In all the
			bead-spring polymer simulations, gamma is used - partially since people
			follow what Kurt and I did. In BD for macro particles, I suspect people
			think of the drag in terms of eta. Maybe easiest to leave a gamma but
			in notes define gamma in terms of eta so people can do the conversion if
			they wish.

			For a LJ fluid, T=epsilson/k_B, eta is of order 1 in dimensionless
			units. (For T=1, pure repulsion, rho=0.85, eta=0.7. Slightly higher
			with r_c=2.5).

			Gary

			Gary - a follow up to the discussion on adding a drag term to LAMMPS
			for implicit solvent effects.

			F_drag = - gamma velocity

			I understand that gamma = 6 pi eta diameter
			where eta = viscosity of the solvent

			which comes from gamma = kT / Dm
			and D = kT / (6 pi eta diameter)

			What would be the typical way a user would want to specify
			this term. To specify the viscosity and then scale gamma for
			different size particles? Do you know what dimensionless viscosity
			is for LJ systems?

			Or just to specify gamma in units of force/velocity for sigma = 1
			particles, and then perhaps scale up gamma for larger/smaller
			particles?

			Any ideas on the normal way to think about the units of this term?

			Steve

			--------------------------------------------
			22 March 2006
			derived classes, rules that are used in LAMMPS
			from talking with Roscoe Bartlett

			assume what is wanted is for child to override some subset of methods
			no matter what calls what, have child's version called if it exists,
			else parent's version

			how to do this in parent class:
			public methods and protected methods
			protected = child can invoke them
			anything that may be overridden by any child is declared a virtual method
			others are not virutal
			no pure virtual = 0 classes if parent is stand-alone
			data (double, int, etc) can be in public or protected section

			how to do this in child (derived class):
			just implement the methods that are new
			can define own variables (but can only refer to them in child class)

			one exception is destructor:
			both child and parent will be called
			is always virtual in parent class
			but parent destructor cannot call child routines (e.g. deallocate)
			at least ones that override parent ones
			this is b/c child may have already gone away
			so parent destructor will always call parent version
			thus need to insure child destructor calls the right thing
			see pair_eam and its FS child for example

			--------------------------------------------
			21 March 2006
			conjugate gradient for eng min

			discussion in Ch 10 of Numerical Recipes is good

			CG requires good line min
			BFGS (quasi-Newton) does not
			stores NxN where N can be small

			CG:
			maintain g = dowhnill gradient (- Grad phi = force vec)
			h = search dir for line min
			basic idea is get a new search dir in outer loop at each iteration
			use h to do linemin to new pt in inner iteration
			linemin does bracketing, then 1s search via Brents method (use deriv info)
			h0 = g0
			doing a linemin gives you new gi+1 at new point
			hi+1 = gi+1 + gamma hi
			gamma = FR or PR forms
			gamma = gi+1 dot gi+1 / gi dot gi for FR

			--------------------------------------------
			14 March 2006
			box vs lattice default

			thought about changing lattice default to box default, but it would
			require lots of people to change their input scripts (i.e. region
			command)

			and code is good now b/c it assumes lattice, so if someone doesn't
			specify a lattice then they get an error

			commands affected:
			velocity (set/ramp), temperature, region, displace_atoms, fix indent

			--------------------------------------------
			10 Jan 2006
			instrumentation to LAMMPS to look at Red Storm performance in detail

			finish.cpp
			add histogram section on pair performance that is in this dir
			to finish.cpp

			atom.cpp, neigh_bond.cpp
			DELTA, BONDDELTA settings are 10000 by default
			get big change on Red Storm if change DELTA

			Makefile.storm (catamount, -D)

			thermo_one.cpp
			change to no output during timestep loop, will get non-discrete
			timings on Red Storm
			just to screen I think, OK to logfile

			timer.cpp
			comment in syncho timings if want

			in.lj
			can change neighbor page and one settings to get different memory
			allocation and performance

			--------------------------------------------
			6 Dec 2005
			use of map array in atom.cpp

			** Routines in atom:

			atom::map_init()
			free old map() array
			allocate map(1:maxtag+1) array
			fill it with -1 on each proc
			called by create_atoms, delete_atoms, read_data, read_restart, replicate,
			velocity (tmp)

			atom::map_clear()
			set map[] to -1 for all my owned and ghost atoms
			called by comm, special (for scratch)

			atom::map_set()
			set map[] for my own + ghost atoms in reverse order
			called by comm, create_atoms, delete_atoms, read_data, read_restart,
			repliacte, special (for scratch), velocity (tmp)

			atom::map_delete()
			free memory for map[]
			called by read_data (tmp), velocity (tmp)

			add int atom::map(i)
			id = global ID
			return -1 or local ID

			** Speed-critical usages of map array:

			comm::exchange()
			at beginning, if (atom->molecular) atom->map_clear()
			at end, if (atom->molecular) atom->map_set()
			dump_bond::count() and pack()
			map lookup of each bond of each atom to check if bonded atom is
			in dump group, if so then global atom ID is output
			fix_shake stores int *map in class so can do lots of lookups
			shake_atom[][] stores global IDs of atoms in SHAKE cluster
			find_clusters() called when fix is invoked
			uses map() to find atoms in bond lists
			pre_neighbor() uses map() to lookup local IDs of those atoms
			also stores atom->map in local int *map
			shake2(), etc uses map() to find atoms in shake_atom[][]
			neigh_bond uses map lookups to find atoms in bond_atom, angle_atom, etc
			and build bond, angle, etc neighbor lists

			** Non-critical usages of map array:

			atom::unpack_vels()
			map lookup of tag of atom in each line to see if I should store vel
			atom::unpack_bonds(), angles, dihedrals, impropers
			map lookup of tag of each atom in bond, etc to see if I should store bond
			atom::memory_usage()
			accounts for memory in map array via maxmap
			create_atoms::create
			at end of create, if atom->molecular: atom->map_init(), atom->map_set();
			delete_atoms::delete
			at end of delete, if atom->molecular
			atom->nghost = 0 - no ghosts, so old ghosts won't be mapped
			atom->map_init(), atom->map_set();
			fix_tmd::readfile()
			uses map lookup to find atoms listed by global ID in TMD input file
			read_data::atoms()
			after atoms are read from data file, if (atom->molecular):
			atom->map_init(), atom->map_set();
			read_data::velocities()
			before vels are read, if (atom->molecular == 0)
			atom->map_init(), atom->map_set();
			after vels are read and assigned:
			if (atom->molecular == 0) atom->map_delete();
			read_restart::read()
			at end of file read and setup, if (atom->molecular)
			atom->map_init(), atom->map_set()
			replicate::
			at end of setup of replicated problem, if (atom->molecular)
			atom->map_init(), atom->map_set()
			set::set()
			uses map lookup to bond_atom, etc to verify that all atoms in a
			bond, angle, etc are in the set group
			special::build()
			uses map lookup to find atoms in bonds and circulate unsatisfied to
			other procs
			special::combine()
			at beginning, atom->map_clear() to allow use of map as scratch space
			at end, atom->map_set() to recreate the map
			velocity::create()
			for loop_flag = ALL, need map to find atoms I own
			at beginning, if (atom->molecular == 0)
			atom->map_init(), atom->map_set()
			at end, if (atom->molecular == 0) atom->map_delete()

			--------------------------------------------
			6 Nov 2005
			Bob Skeel visit and talk

			induced dipoles and fluctuating dipoles are same thing

			he has 2005 JCC paper on fast dipole solver for fluctuating dipole model
			CHARMM has dipole option now

			"shadow" energy is a way to monitor energy conservation very
			accurately - he has recent paper on this

			a symplectic integrator is a phase-space volume-preserving integrator
			energy-conserving is a weaker criterion
			Dan Negrut at U Wisc (formerly Argonne) has worked on this recently
			(for rigid body integrators?)

			--------------------------------------------
			30 Sept 05
			meaning of open source, GPL from a Slashdot discussion

			(Score:5, Insightful) by TrentC (11023) <jelmore.viclink@com> on
			Friday September 30, @03:23AM (#13682722)
			(http://www.crystalwind.org/)

			The GPL, in fact, guarantees that if GPL'd software is used in another
			product, both products then become infected by the GPL and the
			resulting work is then covered by the GPL.

			This is, sadly, a common misunderstanding when it comes to the GPL. By
			using the term "infected", you are either misinformed or attempting to
			misinform; I'll assume the former...

			If you use code licensed by the GPL in your closed-source work and you
			get "caught" distributing it, you have four options:

			1. You can try to obtain an exception to the GPL from the copyright
			holder(s) of the GPLed code for your particular work.

			2. You can change the license of your work to the GPL (or,
			possibly, one of the licenses deemed "GPL-compatible"; IANAL, so
			consult a lawyer first).

			3. You can rewrite the affected portion to remove the GPLed code
			from your work.

			4. You can stop distributing your work so long as it contains the
			GPLed code.

			The copyright holder of the GPLed code can not force you to pick any
			particular one of the options (except, by the definition of the GPL,
			you must do #4 if you can't or won't do #1, #2 or #3). You are the
			copyright holder of your code, and cannot have your license changed
			against your will any more than they can have the license of their
			work changed against their will.

			Jay (=

			--------------------------------------------
			26 Aug 05
			bulliten board ideas

			"Roman Voronov" <burn@ou.edu>
			As far as the message board goes, I like websites that use
			http://www.phpbb.com software(if you click on the Community link you can see
			it at work). However, I have not had any personal experience with running
			it. I would like to see a message board over a mailing list, because it
			would be convenient to share experience with other users, instead of bugging
			you with my emails. Thanks.

			could do sourceforce - Marcus has an account there for Towhee

			--------------------------------------------
			26 Aug 05
			PPPM notes

			when hardwire Gewald and mesh size:
			want mesh spacing (in Angstroms) times Gewald to be ~ 0.25

			error in PPPM part is proportional to Q^2 (h alpha) ^ order
			where Q = size of one Q or size of one Q^2 ??
			h = mesh spacing
			alpha = G_ewald
			order = 5 typically
			so if make charges 10x bigger (due to dipole separation)
			then error goes up by 100, so need to adjust h (in one dimension)
			to compensate

			differences between PPPM and PME and SPME
			PPPM is more accurate than PME for same grid spacing
			SPME reduces the gap, gets closer to PPPM, but not quite as good
			if do 2 FFTs instead of 4 (e.g. in SPME) then lose some accuracy
			but may be more parallel, since less FFTs
			we could do 2 FFTs instead of 4 in PPPM, but with same trade-off
			Gaussian real-space PME may be kind of like multipole
			paper claims it is now as accurate as SPME for 1st time
			but maybe not better yet

			--------------------------------------------
			23 Aug 05
			notes on dipole testing

			erf(x) = 2/sqrt(pi) x as x -> 0

			E-field due to point charge:
			Z = - grad V = - grad (q/r)
			= q/r^3 rij (i.e. it points away from a positive charge
			and falls off as r^2)

			E-field due to point dipole
			Z = - grad V = - grad (1/r^3 p dot rij)
			= 3/r^5 (p dot rij) rij - 1/r^3 p
			so on horizontal axis = -1/r^3 p
			on vertical axis = 2/r^3 p

			** run in 3 modes and compare:

			a) convert each dipole into separated charges with bond
			run with pair lj/cut/coul/long and pppm
			b) convert each dipole into separated charges with bond
			run with pair dipole/long and pppm/dipole
			c) treat each dipole as true point dipole
			run with pair dipole/long and pppm/dipole

			compare energy, force, torque on each particle (or pair of particles)

			** energy correction factor for one dipole in a periodic box

			e_ij = qi qj / r
			f = weighting factor on a bond = 0 to 1
			when no bond: E_nb = e_ij erfc + E_long
			with bond: E_b = E_nb - e_ij (1 - f)
			when f = 0: E_b = e_ij (erfc - 1) + E_long
			true dipole: E_d = 0 (short range) + E_long (with separated charges)

			our dipole is like 2 charges with f = 0 bond (full exclusion of charge-charge)
			so the E_long for dipole we compute in pppm_dipole.cpp
			needs to reproduce E_b formula when f = 0
			means for each dipole, should add in e_ij (erfc - 1)
			to match short range contribution of 2 point charges

			this is done in init() of pppm_dipole.cpp to compute eng_correct

			** other energy and virial correction factors

			for energy in pppm_dipole:

			** one dipole in 20^3 box, sep = 0.1,
			PPPM accuracy = 0.0001, mesh = 64^3, gewald = 1.0

			a) Etot = -0.00046063816, Ecoul = 112.463
			b) Ecoul = 112.46297, E_long = -112.46343, Etot = -0.00046063816
			c) Ecoul = 0.0, E_long = Etot = -0.00046063816
			E_pppm = -112.463, E_correct = 112.463

			so good agreement

			** two dipoles in 20^3 box, sep = 0.1, placed at -1 and 1 in x
			PPPM accuracy = 0.0001, mesh = 64^3, gewald = 1.0, RMS = 0.00546

			a) Etot = -0.25182084, Ecoul = 224.872
			Fx = -2.33651736274 + 2.71299369315 = 0.376476330410
			b) Ecoul = 224.87238, E_long = -225.1242, Etot = -0.25181978
			Fx = -2.33652035228 + 2.71299585767 = 0.376475505390
			c) Ecoul = -0.052836928, E_long = -0.19826378, Etot = -0.25110071
			Epppm = -225.124, E_correct = 224.926
			Fx = 0.374021405918

			good agreement, better in eng than force
			could disagreement be in short-range, not long-range ??

			** two dipoles in 20^3 box, sep = 0.02, placed at -1 and 1 in x
			PPPM accuracy = 0.0001, mesh = 128^3, gewald = 2.0, RMS = 0.00413

			a) Etot = -0.21246041, Ecoul = 11277.8
			Fx = -13.1373720414 + 13.5124040759 = 0.37503203449999845
			b) Ecoul = 11277.811, E_long = -11278.024, Etot = -0.21246041
			Fx = -13.1380970108 + 13.5131290459 = 0.37503203510000027
			c) Ecoul = -2.162565e-06, E_long = -0.21245823, Etot = -0.21246039
			Fx = 0.375031760716

			even better agreement, but could be b/c of higher gewald
			pushing more computation into long-range

			** one dipole at -1, one pt charge at +1
			dipole on top of point charge
			two dipoles with 2 point charges

			good agreement in all cases for energy and forcesbetween 3
			virial close but small so not perfect agreement

			need to check on virial as run bigger systems

			bugs I found: 1/2 factor in separation in pppm_dipole
			didn't have q and dipole included in rho_map and unmap
			needed E correction due to paired-charge term no longer computed in
			pairwise

			--------------------------------------------
			21 Jun 05
			Seel notes on LJ + EAM via hybrid vs modifying EAM to do LJ:

			You'll be happy to know that the LJ + EAM hybrid capability gives the
			exact same energy/pressures as doing the same with a SETFL, but the
			hybrid method is 30% faster. In order to get exactly the same answers
			(mostly because of the funniness of putting a LJ into a SETFL) I had
			to use shifted and smoothed (ie., derivative=0 at rcut) LJ potentials
			for both the hybrid and SETFL runs, which necessitated using
			"pair_style table". I've attached all the input scripts, data files,
			and output dumps for your enjoyment.

			--------------------------------------------
			21 Jun 05
			notes on NVT and NPT

			FS = Frenkel Smit book on Understanding Molecular Simulation
			Mel = Melchionna 1993 paper on NPT w/ center-of-mass
			Mar = Martyana 1995 paper on NVT/NPT in exp() notation for hierarchical steps

			** Nose/Hoover NVT (FS eq 6.1.24-6.1.27, Mel eq 7 w/out pressure part)

			dx/dt = v
			dv/dt = F/m - eta v
			d(eta)/dt = v_T^2 (T/T_ext - 1)

			** Verlet form of N/H NVT (straightforward differencing)

			eta' = eta + dt/2 v_T^2 (T/T_ext - 1)
			v' = v + dt/2 (F/m - eta' v)
			x_new = x + dt v'
			F_new = force(x_new)
			v_new = v' + dt/2 (F_new/m - eta' v')
			T_new = temperature(v_new)
			eta_new = eta' + dt/2 v_T^2 (T_new/T_ext - 1)

			** NVT with exp() operators (Mar eqs 26,35 and Appendix A)

			no chains, set NNOS = w = nresn = nyosh = 1, Nf = 3N deg of freedom

			is just 2 scalings of v on either side of NVE velocity Verlet

			v = v * scale
			v' = v + dt/2m F
			x_new = x + dt v
			F_new = force(x_new)
			v_new = v' + dt/2m F_new
			v_new = v_new * scale

			scale computation:

			KE = sum (m v^2)
			GLOGS = (KE - 3NkT_ext) / Q
			VLOGS += dt/4 GLOGS
			scale = exp(-dt/2 VLOGS) (roughly 1 - dt/2 VLOGS)

			** LAMMPS NVT (fix_nvt.cpp)

			f_eta = v_T^2 (T/T_ext - 1)
			eta_dot' = eta_dot + dt/2 f_eta
			factor = exp(-dt/2 eta_dot') (roughly 1 - dt/2 eta_dot')
			v' = v * factor + dt/2m F (roughly v' += dt/2m F - dt/2 eta_dot')
			x_new = x + dt v'

			F_new = force(x_new)

			v_new = (v' + dt/2m F_new) * factor (factor applied to both terms)
			T_new = temperature(v_new)
			f_eta = v_T^2 (T_new/T_ext - 1)
			eta_dot_new = eta_dot' + dt/2 f_eta

			** Comparison of all NVT options

			LAMMPS NVT is similar to Verlet form of N/H
			if exp(-dt/2 eta_dot') = 1 - dt/2 eta'
			which it should be for small eta_dot
			will eta_dot ever accumulate (after many timesteps) to large value
			only if T_actual stays larger (or smaller) than T_desired for a long time
			LAMMPS applies old factor to (v' + F)
			maybe that is like applying eta_new in Verlet difference ??
			Verlet has ambiguity in that eta_new requires T_new which requires
			v_new which requires eta_new - so not sure which to use when
			LAMMPS NVT is similar to Martyana NVT with exp() operators
			if VLOGS = eta_dot
			so GLOGS/2 = f_eta
			and KE - 3NkT_ext / 2Q = v_T^2 (T/T_ext - 1)
			since sum (m v^2) = KE = 3NkT
			thus v_T^2 = 3/2 NkT_ext / Q
			Martyana (p 1121) says Q = 3NkT / w^2 where w is freq of T oscillations
			so v_T^2 = w^2 / 2
			not sure if w is frequency or radians/sec
			but if former, the 2 methods match within a factor of sqrt(2)
			one difference: Martyana computes KE (T_new) in 2nd scale based on current v,
			which means it has been updated by F_new but not 2nd eta, that
			quantity is never available in LAMMPS b/c v is updated all at once

			** Nose/Hoover NPT (Mel eq 7), in center-of-mass form

			dx/dt = v + omega (x - X0) (X0 is center-of-mass)
			dv/dt = F/m - (eta + omega) v
			d(eta)/dt = v_T^2 (T/T_ext - 1)
			d(omega)/dt = v_P^2 (Vol/NkT_ext) (P - P_ext) (Vol/NkT is like 1/P_ext)
			d(Vol)/dt = 3 Vol omega

			** Verlet form of N/H NPT (straightforward differencing)

			eta' = eta + dt/2 v_T^2 (T/T_ext - 1)
			omega' = omega + dt/2 v_P^2 (Vol/NkT_ext) (P - P_ext)
			v' = v + dt/2 (F/m - (eta' + omega') v)
			Vol' = Vol + dt/2 (3 Vol omega')
			x_new = x + dt v'

			F_new = force(x_new)

			Vol_new = Vol' + dt/2 (3 Vol' omega')
			v_new = v' + dt/2 (F_new/m - (eta' + omega') v')
			T_new = temperature(v_new)
			eta_new = eta' + dt/2 v_T^2 (T_new/T_ext - 1)
			P_new = pressure(x_new,Vol_new,virial_new,T_new)
			omega_new = omega' + dt/2 v_P^2 (Vol_new/NkT_ext) (P_new - P_ext)

			** NPT with exp() operators (Mar eq 39 and Appendix B)

			no chains, set NNOS = w = nresn = nyosh = 1, Nf = 3N deg of freedom

			is just 2 scalings of v on either side of NVE velocity Verlet
			plus box rescaling of coords and volume

			v = v * scale
			v' = v + dt/2m F
			AA2 = exp(dt VLOGV)
			x_new = x * AA2 + dt v (actually modified more for box velocity)
			F_new = force(x_new)
			v_new = v' + dt/2m F_new
			v_new = v_new * scale

			scale computation:

			KE = sum (m v^2)
			GLOGS = (KE + VMASSVLOGVVLOGV - (3N+1)kT_ext) / Q
			GLOGV = [(1 + 1/N) * KE + 3 (P_int - P_ext) * Vol ] / VMASS
			VLOGS += dt/4 GLOGS
			AA = exp(-dt/8)
			VLOGV = VLOGV AA AA + dt/4 GLOGV AA
			AA = exp(-dt/2 (VLOGS + (1 + 1/N) VLOGV))
			scale = AA
			KE = KE AA AA
			GLOGV = [(1 + 1/N) * KE + 3 (P_int - P_ext) * Vol ] / VMASS
			AA = exp(-dt/8 VLOGS)
			VLOGV = VLOGV AA AA + dt/4 DLOGV AA

			** LAMMPS NPT (fix_npt.cpp)

			f_eta = v_T^2 (T/T_ext - 1)
			eta_dot' = eta_dot + dt/2 f_eta
			f_omega = v_P^2 (Vol/NkT_ext) (P - P_ext)
			omega_dot' = omega_dot + dt/2 f_omega
			factor = exp(-dt/2 (eta_dot' + omega_dot'))
			dilation = exp(dt/2 omega_dot')
			v' = v * factor + dt/2m F
			xprd' = xprd * dilation (box dilation)
			x_new = x + dt v'

			F_new = force(x_new)

			xprd_new = xprd' * dilation (box dilation)
			v_new = (v' + dt/2m F_new) * factor (factor applied to both terms)
			T_new = temperature(v_new)
			f_eta = v_t^2 (T_new/T_ext - 1)
			eta_dot_new = eta_dot' + dt/2 f_eta
			P_new = pressure(x_new,xprd_new,virial_new,T_new)
			f_omega = v_P^2 (Vol_new/NkT_ext) (P_new - P_ext)
			omega_dot_new = omega_dot' + dt/2 f_omega

			--------------------------------------------
			28 Apr 05
			new NPT ideas for better pressure equilibration

			see NOTES/npt_ideas
			plots from Jon Z of problems with NPT Nose/Hoover for metals
			2 fix_npt_pr files he started to allow Parinello/Rahman algorithm instead

			--------------------------------------------
			14 Mar 2005
			how EAM files are read/interpolated in original LAMMPS 2004 version

			simulation is either all funcl, or all setfl (single file, all types mapped)

			** funcfl (one element):

			read_funcl: once for each pair_coeff command
			mass, nrho, drho, nr, dr, cut
			arrays: frho(nrho), zr(nr), rhor(nr)

			if i = j (in pair_coeff)
			set global mass
			store_funcl:
			stores read-in values in [type][nr or nrho] arrays
			so this is done just for each type

			convert_funcfl (called once by init for all ntypes):
			cutforce = max of all cuts
			ditto for dr,drho
			set global nr,nrho based on dr,drho and longest potential range
			interp frho[itype][nrho]
			interp rhor[itype][nr]
			interp zrtmp[itype][nr] from zr
			set z2r[type][jtype] = 27.20.529 zrtmp[i] * zrtmp[j]

			spline coeff interp from frho,rhor,z2r

			set cutsq[i][j] to cutforce^2
			cutforcesq = cutforce^2, used in compute() for all pairs

			** setfl (multi-element):

			read_setfl: one time
			mass[itype], nrho, nr, drho, dr, cut
			arrays: frho[itype], rhor[itype], z2r[itype][jtype] for itype <= jtype
			copy z2r values for itype > jtype

			create map[itype] = arg in pair_coeff command
			set global mass of all types

			global nr,nrho,dr,drho,cutforce set from setfl file values

			copy frho,rhor,z2r file values into arrays used by interpolate()

			spline coeff interp from frho,rhor,z2r

			set cutsq[i][j] to cutforce^2
			cutforcesq = cutforce^2, used in compute() for all pairs

			--------------------------------------------
			25 Jan 2005
			how to scan directory for filenames

			works on Linux, Tflop, should soon work on other machines

			#include <stdio.h>
			#include <sys/types.h>
			#include <dirent.h>

			int main (void)
			{
			DIR *dp;
			struct dirent *ep;

			dp = opendir ("./");
			if (dp != NULL)
			{
			while (ep = readdir (dp))
			puts (ep->d_name);
			(void) closedir (dp);
			}
			else
			perror ("Couldn't open the directory");

			return 0;
			}

			--------------------------------------------
			25 Jan 2005
			PPPM and PME discussion with Paul

			used to be that Ewald used beta (Ewald G-vector) and PPPM used G/sqrt(2)
			see Roy paper on PPPM - uses G/sqrt(2)

			Paul changed it so that both Ewald and PPPM just use G
			consistent with PME papers of Darden, etc
			did this by changing how G is set in PPPM setup

			PME formulations can be done to do 2 FFTs instead of 4 like we do in PPPM
			done by analytic form of derivatives
			Darden dipole paper has some discussion of PME variants like this
			Paul has paper that carefully compared PME/PPPM and these variants
			for charge system (no dipoles)
			concluded PPPM is better than PME
			smoothing can be done for PME or PPPM (in different way)
			trade-off in accuracy vs FFTs
			4 FFTs is more accurate so fair comparison is against 2 FFTs of larger grid
			in parallel, the FFTs hurt, so may be better to do 2 FFT version

			--------------------------------------------
			18 Jan 2005
			Pair routine and bond routine forces

			r = scalar distance between atoms i and j
			rij = vector = R_i - R_j, delx = x_i - x_j = vector from j to i
			E(r) is given for pair, bond, etc
			F = -Grad_i E(r) = - dE/dr Grad r
			Grad r = dr/dx i + dr/dy j + ...
			dr/dx = 1/r (x_i - x_j)
			so Grad r = 1/r rij
			so F = -1/r dE/dr rij

			in pair routines:
			forcelj variable (lj/cut) = - r dE/dr
			fforce variable (all pair routines) = - 1/r dE/dr
			delx = x_i - x_j
			Fij_x = x-force on atom i due to atom j = fforce * delx

			in bond routines:
			fforce variable = - 1/r dE/dr

			--------------------------------------------
			17 Jan 2005
			Releasing a new version of LAMMPS

			do these steps in order

			lammps/src:
			pur (to get rid of old files)
			make package-check to insure all src files are current with
			master copies in package sub-dirs
			take out any changes/files/packages that are not ready for release
			include style files and package/Install.csh in this check

			tests:
			don't have regression set yet - so skip for now
			idea: run all tests and compare to old files to see if they work

			lammps/tools:
			is restart2data up-to-date with all potentials and atom style in src
			is replicate.c up-to-date with any new atom styles

			lammps/doc:
			pur (to get rid of old files)
			txt2html *.txt
			Section_intro.txt
			update feature sub-section (see WWW page version)
			Section_history.txt
			update future plans sub-section (see WWW page future list)
			Setion_errors.txt
			run error.py to insure error list and source code agree
			make sure all packages are included in src for this operation
			browse Section_errors.html to insure haven't left off a tag
			make sure all Eqs are up-to-date by running Script (update Script as well)
			if necessary, re-import to JPG, clip, at 33% from xpdf display

			lammpsdist.sh 17Jan05 '17 Jan 2005' (one year is 2-digit, other is 4)
			if are any new packages, examples, bench, tools, doc, lib files, etc then
			edit bin/lammpsdist.sh to include them
			include ALL packages in src so will be able to run benchmarks & examples
			run lammpsdist.sh to create a new tar file temporarily
			this will change date in all src files before running bench & examples
			will need to re-create this tar file after run bench & examples

			bench:
			create a new bench/Log dir with the new date
			update Script.compare, Script.pack, Script.unpack
			make new executables for linux and liberty, etc
			run 5 benchmarks on 1 and P procs, using new executables
			run Script.pack to pack benchmark files for any machine
			edit and use Script.linux.1, nqs.liberty.8, etc
			put log file results in new Log_date dir
			run Script.unpack to put those log files back into Lj, Eam, Chute, etc
			run Script.compare to see if answers/timing changed
			log.date files in Lj, Chute, Chain, Eam, Rhodo are ones copied to
			new dist, so make sure new dated copies are in those 5 dirs

			examples:
			update Script.compare, Script.pack, Script.unpack
			run Script.pack to pack example files into HOLD, xfer to whatever machine
			re-run all examples with new version via Script.linux.1 and Script.liberty.4
			set date and executable in scripts first
			need to run liberty one as interactive batch else change dir seems to fail
			copy log* and dump* files back from parallel machines into HOLD
			run Script.unpack to put log* and dump* files back into example dirs
			run Script.compare to see if answers/timing changed
			viz any movies from new dump files if worried answers not the same
			gzip all the dump* files
			all log.date files will get put in distribution
			clean up and make sure they are current with new version

			lammpsdist.sh 17Jan05 '17 Jan 2005' (one year is 2-digit, other is 4)
			final build of tar file with all examples and bench runs
			first whack previous tar file and dir in lammps/versions
			make no-PACKAGE for any packages not to be included by default
			cp tar file into tmp dir, unpack and build LAMMPS to make sure it works
			look thru unpacked dir to insure everything is there
			new packages, new examples, new versions of README files, etc
			check new Manual.pdf
			add list.17Jan05 as empty file to lammps/patches
			make sure package dirs are listed correctly (including new ones) in
			bin/patch.py

			2001:
			make new release if necessary
			save old tar file in versions/2001 with date
			create tar file in versions/2001
			put date of new release in download.txt and size of new 2001 tar file

			text/www/
			index.txt
			set new date for current LAMMPS version
			create a "New.gif" bullet for the new release
			change hi-lite to a new picture or movie, test all links
			abstract and pic links need to be different
			store line in old hi-lites at bottom
			bug.txt
			set new date for current LAMMPS version
			collapse all old patches into a file like bug.1Sept04
			make an entry on release date with list of all major changes
			put in line count for new version
			download.txt
			set new date for current LAMMPS version (2)
			set size of new tar file
			save WWW:download.log file on my machine
			update count of downloads
			count.py download.log LMP 600 0 (10 min delay)
			count.py download.log.2004-6 LMP -1 0 (infinite delay)
			also add in SF downloads to count
			also do ALL to get non-LAMMPS downloads for my WWW download page
			features.txt - update it (see Section_intro of doc)
			future.txt - update it (see Section_history of doc)
			txt2html *.txt in both text/www and text/www/lammps

			actual update:
			www.py -dist lammps 17Jan05
			see www.py header for full info
			puts all doc files in text/www and onto WWW server
			puts tar file and new soft-link on WWW server
			on WWW in tars: cp lammps-17Jan05.tar.gz lammps-upgrade.tar.gz
			text/www
			www.py -p download.txt download.html lammps.txt lammps.html lammps_new.*
			text/www/lammps
			www.py -p .txt .html
			put new 2001 tar file onto barky if changed
			try to do a download from WWW site myself
			put new tar file on SourceForge
			click on "file releases" under admin
			click on "add release" at bottom of page
			new release name = lammps-3Jun05
			no release notes or change log, click submit
			goto Step 2
			upload tar file via anon ftp and select it to add the file to release
			ftp to upload.sourceforge.net, cd to incoming, upload the file
			click on refresh to see new file, select it, click on add file
			should see File Updated
			goto Step 3
			processor release date = Any
			file type update = Source .gz
			click on update/refresh
			can then test a download
			send email to lammps-users mail list if desired via Gmail

			--------------------------------------------
			3 Jan 2005
			fix rigid integration and accuracy vs SHAKE

			see NOTES/fix_rigid
			graphs of temperature, energy, pressure for 2 runs of same system
			a box of a few 100 waters, held rigid by SHAKE or fix rigid
			pressure does not agree on 1st timestep

			bad to integrate with just normal NVE w/ velocity Verlet
			energy is not conserved
			issue of when to renormalize Q
			issues of dq/dt = 1/2 q w_body
			w_body depends on q, so aren't being consistent about where
			in timestep each quantity is when do update
			think it becomes 1st order with these errors, hence very bad energy conserve
			options: Richardson extrapolation or Nic Martys solution
			we used Richardson because it was cheap and as good as an iterative
			solver in Nic Martys paper
			is also a question of when to stop iterations

			put code for other fix rigid integrators we didn't end up using
			in NOTES/fix_rigid
			vanilla = old bad one
			iterate = Nic Martys style one

			Paul's tests that justify use of Richardson

			Here are the energy conservation results.

			Method log(ave abs of energy fluctuation)
			shake -3.913702852
			vanilla -0.625317489
			rich -3.538767441
			iterate.2 -3.48685072
			iterate.3 -3.518183999
			iterate.4 -3.569168296
			iterate.5 -3.587866192
			iterate.10 -3.573231918

			Timings were essentially the same. This was for a 10 ps run of 216
			waters. My conclusions:
			1) shake gives the best energy conservation
			2) all the others except vanilla give good energy conservation
			3) more than 4 iterations doesn't substantially improve E conservation

			So, I vote for releasing it with either rich or 4 or 5 on the iteration
			method.

			--------------------------------------------
			3 Jan 2005
			Paul's graphs for bit-mapped table performance for long-range Coulombics

			see NOTES/tables_coulomb
			graphs of CPU performance and RMS error bounds

			--------------------------------------------
			22 Dec 2004
			Dipoles

			ways that dipole length is used:

			atom_dipole::init() - normalizes each atom's dipole
			dump - just dumps it
			integrator nve - renormalizes each atom to dipole length
			set_type - creates a dipole of proper length
			initial data file - just reads it
			create atoms - sets dipole to zero

			--------------------------------------------
			14 Dec 2004
			Force tabling for coul long

			code for how I did lookup, linear, spline before Paul chose just
			linear with bitmap is in NOTES/table_coulomb

			--------------------------------------------
			7 Dec 2004
			Full neighbor lists - and how to use them

			added neigh_full routines to construct neighbor lists where
			all the neighbors of each atom are listed

			changes to neighbor class:
			init:
			full flag is 0 if only one regular list is maintained
			full flag is 1 if only one full list is maintained
			full flag is 2 if an extra full list is maintained
			looks at fixes, dumps to decide what to set full flag to
			2 versions of neigh_full routines: nsq and bin
			neigh_full works using regular neigh list variables (pages, firstneigh, etc)
			build_full calls nsq or bin version of neigh_full, but swaps
			regular vars with full vars (pages_full, etc) so neigh_full constructs
			list in extra memory
			can construct full list in regular memory by setting pair_build fn ptr
			to neigh_full nsq or bin
			setup_bins constructs a stencil for full neighbor list if full
			rest of neighbor routines manages memory for 2 sets of list if full flag = 1

			used by Si potentials (3-body)
			only 1 neighbor list: full
			call directly from neighbor->build() directly by setting ptr
			used by Koen in his fix to compute extra forces due to atom neighborhood
			he will want updated list every timestep
			rest of LAMMPS needs normal list
			neighbor::init() looks for named fix (koen) to set full flag = 2 (extra list)
			fix calls neighbor::build_full() when neighbor->ago is 0
			fix uses firstneigh_full and numneigh_full to access full list
			used by dump of centro symmetry (CENTRO fix)
			fix needs full list every 100 steps (or thereabouts)
			rest of LAMMPS needs normal list
			neighbor::init() sets full flag = 2 if a fix centro is defined
			fix calls neighbor::build_full() each time it is invoked at end-of-step
			fix uses firstneigh_full and numneigh_full to access full list
			will be inefficient if 2 or more dumps use CENTRO
			will be inefficient if dump is called (nearly) every step

			--------------------------------------------
			26 Oct 2004
			Changes vs F90 LAMMPS

			in F90, force units in dump are not good
			in F90, -3 dof in temp is not done in velocity creation, nor is SHAKE dof
			do it right in 2003
			in SHAKE, the shake_update in src F90 version uses full dt for
			unconstrained update even when dtsq_factor is 1/2 - this seems
			incorrect
			use of PI in code is inexact in F90 LAMMPS, look for 3.14159 thruout code
			shouldn't need extra param in Temp creation
			# of degrees of freedom to subtract
			Langevin should not be pegged to Marsaglia RNG
			way fix_Langevin calls initial seed with proc ID
			langevin is different in LAMMPS 2001 when done with fix via temp rescale
			in former case, no call is done from start routine (itime = 0)
			not sure which is correct for restarting - in new LAMMPS there is
			only one kind of langevin fix, so need to make pre() correct for both
			ways - with call in pre(), matches temp rescale in LAMMPS 2001,
			without call in pre(), matches fix langevin in LAMMPS 2001
			think about which way is correct!

			--------------------------------------------
			25 Oct 2004
			User question about LAMMPS

			Here is the kind of "help me" email you get when you put
			your code on the WWW for anyone to grab ...

			Maybe I'll tell him to send mail to Microsoft asking why
			there isn't a "run" command in Linux ...

			I already built LAMMPS by=20
			"make serial"
			command and then I got lmp_serial fiel in src directory.
			Then I copy it to example/melt directory and command=20

			"run lmp_serial < in.melt".

			But I only got a message
			"bash: run: command not found".

			--------------------------------------------
			15 Oct 2004
			What GPL license means for another code:

			HTMLDOC is provided under the GNU General Public License ("GPL") with
			a license exception for the OpenSSL toolkit. A copy of the exception
			and license follows this introduction.

			For those not familiar with the GNU GPL, the license basically allows you to:

			* Use the HTMLDOC software at no charge.
			* Distribute verbatim copies of the software in source or binary form.
			* Sell verbatim copies of the software for a media fee,
			or sell support for the software.
			* Distribute or sell your own modified version of HTMLDOC so long as
			the source code is made available under the GPL.

			What this license does not allow you to do is make changes or add
			features to HTMLDOC and then sell a binary distribution without source
			code. You must provide source for any changes or additions to the
			software, and all code must be provided under the GPL.

			--------------------------------------------
			10 Sept 2004
			EAM in LAMMPS vs Warp vs ParaDyn

			why is ParaDyn so slow

			Here are some quick timings (CPU secs for 100 steps) I did for this
			32K atom problem on Tflops (ross is down, I don't think it would be
			much different, just faster):

			Procs ParaDyn LAMMPS Warp
			1 618 176 155
			8 73.6 23.9 20.3
			64 10.4 3.60 3.03

			All the codes scale OK. ParaDyn probably gets a bit of a boost from
			cache usage as the speed-up from 1->8 procs is unrealistic. Warp is a
			little faster than LAMMPS, b/c it is a much more streamlined code.

			The real issue is how much slower ParaDyn is on one-proc. This goes
			back to when Warp was written - both it and Pdyn are Fortran (f77 and
			f90). I thought I just took ParaDyn code for force evaluation and
			stuck it into Warp, but I went back and looked at my Warp notes and
			benchmarks. It was 3-4x faster than ParaDyn on all machines from the
			start. Don't remember why. I must have done something faster than
			ParaDyn with how neighbor lists are used and forces are computed in a
			pairwise sense. That carried over into LAMMPS as I lifted code from
			Warp for LAMMPS (converting to C++, which doesn't change performance
			much).

			--------------------------------------------
			20 Aug 2004
			Tool for auto-generate of a distribution

			** tool that you give distribution date to

			set version_str in universe.cpp
			set new date in Manual.txt
			set date on download page
			set date on bug page
			anyplace else in DOC or lammps WWW site that date is?
			run txt2html on doc dir
			copy Doc files to lammps/doc
			create a new Makefile.explicit

			** create Distribution dir and tar it:

			lammps_4Sept2004

			README - look in doc - you should have N dirs and 2 files - point to LWS - go to Manual.html
			go to Sec_start.html, Sec_intro to learn about LAMMPS, my email
			LICENSE
			potentials
			README - DYNAMO format - no DOC - see pair_style EAM
			cuu3
			bench
			README - how to run them, send them to me - see Benchmaarg Page of LWS
			in.eam
			data.chain
			in.chain
			in.chain.P
			log.chain.linux.1
			log.chain.tflop.8
			all 5
			doc
			Eqs - just gif
			Manula.html
			Sec.html
			command.html
			all *.txt
			src
			MAKE/Makefile.*
			Makefile
			Makefile.explicit
			*.cpp
			*.h
			STUBS dir
			examples
			README - see Section_examples - how to run each of 8 (plus 2)
			flow
			in.flow
			log.flow.linux.1
			log.flow.tflop.8
			friction
			all 8
			tools
			Makefile
			top-level tools
			dirs of package tools

			** set of WWW pages

			top-level
			main.html
			*.html
			lammps.html
			download.html
			download.pl
			abstracts
			papers
			pictures
			movies ?
			doc (for all but LAMMPS)
			tars (tarballs) - for each download
			*.txt for each html
			lammps
			FAQ.html
			*.html
			*.txt
			abstracts
			movies
			pictures
			doc - identical to tarball
			inputs
			for bench
			for examples

			--------------------------------------------
			14 Aug 2004
			Aidan comments on NPT in LAMMPS

			I went through the derivation of the Nose-Hoover style barostat
			and obtained a slightly different version than what is in
			LAMMPS. When I tried both versions in GRASP, the LAMMPS version
			showed poorer conservation of the psuedo-Hamiltonian.
			I pointed this out to Mark Stevens, but he did not pursue it.
			I really don't know if it is a big deal, since there are other issues
			with those equations of motion anyway. But I think it
			would be good to report conserved quantities in LAMMPS,
			as a sanity check for users running extended-system simulations.

			--------------------------------------------
			5 Aug 2004
			Mark Stevens interest in large LAMMPS simulations on Red Storm

			Note: the adhesives sim would require bond breaking and formation

			I have some interest in running some large adhesive simulation.
			With Tom Woolf, I am interested in doing some large scale simulations of
			membrane fusion.

			>Mark & Gary -
			>
			>2 things I've heard about the new Red Storm machine that may be of
			>interest to you.
			>
			>(1) 1/4 of it (10 Tflops) should be up and running by Oct. The full
			>machine not until next spring or summer.
			>
			>(2) There will be a bias in the policies towards large jobs (1000s of
			>procs). Since the individual procs will also be fast, this translates
			>into large LAMMPS runs. So you might want to be thinking about big
			>problems that would be cool to try out before the machine gets
			>overloaded. Also, the new LAMMPS has the option to grab 1000 procs
			>and run 100 10-proc jobs under the covers, but that probably isn't what
			>they have in mind ...
			>
			>Steve
			>
			>

			--------------------------------------------
			3 Aug 2004
			broken LADERA link

			Under Publications, the link to LADERA is broken.
			Dynamics of Exchange at Gas-Zeolite Interfaces I: Pure Component Butane
			and Isobutane, M. Chandross, E. B. Webb III, G. S. Grest, M. G. Martin,
			A. Thompson, and M. Roth, J Phys Chem B 105, 5700 (2001). (this paper
			used LADERA <http://www.cs.sandia.gov/~gsheffe/gcmd/gcmd.html> , a MC/MD
			code which uses LAMMPS for its MD portion).

			--------------------------------------------
			21 Jun 2004
			Paul's response to rRESPA

			Looks great!

			I don't see a problem with having SHAKE constrained at every level of
			respa, because typically, our innermost time step would be 2.0 fs if
			we're running SHAKE, and that's the same as the non-respa time step that
			I'm currently using. It won't be any slower. All it will do is allow us
			to take a bigger time step on the slowly changing long-range forces,
			which will allow the simulation to proceed faster.

			The fix rigid we'll probably need to change so that it don't change
			positions, just velocities in the spirit of that paper I sent. We can
			think of it in terms of an impulse "kick" that pushes the bodies by
			altering the velocities of the individual atoms that make up the body.

			I like the idea of the outermost timestep being the one that the user
			specifies and having the neighborhood lists trigger off of it. The inner
			neighbor lists could either be updated every time the outer forces are
			computed, or else have a trigger of their own. But for efficiency, we'll
			probably want to update the inner lists only when the outer forces are
			computed.

			If the langevin forces are inexpensive and need to be updated
			frequently, the inner only may be the right place for them? In fact, any
			of these that are cheap could be done on the innermost level only,
			without substantial slowdown.

			I'm not sure on the NPT and NVT ensembles. It may be adequate to do
			these on the innermost level only.

			Paul

			--------------------------------------------
			21 Jun 2004
			rRESPA notes

			** Input specification

			run_style verlet # the usual velocity-Verlet
			run_style respa N 2 2 2 2 angle 2 pair 4 kspace 5

			N = # of levels
			2 2 2 2 = loop factor between levels, specify N-1 values
			rest or args are optional - assign a force to a level (1 to N)

			there are 6 kinds of forces to assign to levels
			bond, angle, dihedral, improper, pair, kspace
			(could add more like pair/inner, pair/outer, etc)
			by default bond is assigned to 1st level (inner-most)
			by default pair is assigned to Nth level (outer-most)
			be default angle is assigned to bond's level
			be default dihedral is assigned to angle's level
			be default improper is assigned to dihedral's level
			be default kspace is assigned to pair's level

			timestep 10.0 # dt for outermost level
			run 1000 # 1000 steps of outermost level

			** Fixes in rRESPA

			each level has 3 places to call fixes:

			INITIAL_INTEGRATE (typically where v,x are updated)
			POST_FORCE (right after forces are computed)
			FINAL_INTEGRATE (typically where v is updated)

			additionally, the entire hierarchy has these calls:

			END_OF_STEP (at outer-most level)
			PRE_NEIGHBOR and PRE_EXCHANGE (before neighbor lists are rebuilt)

			** Where to apply fixes in rRESPA

			*addforce = POST_FORCE, outer only
			*aveforce = POST_FORCE, every level
			do ave at every level, only add in extra force at outermost
			*com = END_OF_STEP
			*drag = POST_FORCE, outer only
			*enforce2d = POST_FORCE, every level
			*freeze = none (no respa with granular)
			*gravity = none (no respa with granular)
			*gran/diag = none (no respa with granular)
			*indent = POST_FORCE, outer only
			*insert = none (no respa with granular)
			*langevin = POST_FORCE, outer only
			*lineforce = POST_FORCE, every level
			*msd = END_OF_STEP
			*npt = INITIAL & FINAL INTEGRATE, every level
			*nvt = INITIAL & FINAL INTEGRATE, every level
			*nve = INITIAL & FINAL INTEGRATE, every level
			*nve/granular = none (no respa with granular)
			*planeforce = POST_FORCE, every level
			*setforce = POST_FORCE, every level,
			set f = 0 on all inner, set to full value on outermost only
			*shake = PRE_NEIGHBOR, whenever pair-wise neighbor lists re-built
			POST_FORCE at innermost level
			*spring = POST_FORCE, outer only
			*temp/rescale = END_OF_STEP
			*viscous = POST_FORCE, outer only
			*vol/rescale = END_OF_STEP
			*wall/gran = none (no respa with granular)
			*wall/lj93 = POST_FORCE, outer only
			*wiggle = POST_FORCE, outer only

			--------------------------------------------
			18 Jun 2004
			Debye-Huckle notes

			E(r) = q q / r * exp(-Kr) where K = kappa = user input coeff
			used in place of standard Coulomics
			usual leading dielectric coeff
			used with typical additional LJ term

			F = - Grad (E)
			in pair routines: fforce = -1/r * dE/dr
			so that Fxi = x force on atom i = delx*fforce where delx = xi - xj
			gives correct sign of forces for 2 +q charges separated in x,
			with x1 to the right of x2, i.e. atom 1 has +Fx, atom 2 has -Fx
			since dE/dr = -q q / r^2

			E(r) = q q / r * exp(-Kr)
			dE/dr = q q [ -K exp(-Kr) / r + exp(-Kr) (-1/r^2) ]
			= - q q K / r exp(-Kr) - q q / r^2 exp(-Kr)
			= - q q exp(-Kr) [ K/r + 1/r^2 ]
			fforce = q q / r^2 exp(-Kr) [K + 1/r]

			--------------------------------------------
			8 Jun 2004
			Langevin units notes

			force = gamma_1 velocity + gamma_2

			gamma_1 = units of mass / time
			force->ftm2v converts ft/m to v
			thus f_1 = gamma_1 / ftm2v * velocity
			divide gamma_1 in code by ftm2v

			gamma_2 = units of sqrt(m) sqrt(E/t^2)
			force->mvv2e converts mv^2 to eng
			thus gamma_2 = mv/t sqrt(1/mvv2e)
			thus f_2 = gamma_2 / ftm2v / sqrt(mvv2e)
			divide gamma_2 in code by ftm2v*sqrt(mvv2e)

			--------------------------------------------
			8 Jun 2004
			Units notes

			for lj, all conversions = 1
			for real & metal, convert both to cgs and set equal, solve for c

			(0) Boltzmann factor

			should be set to eng/temperature for that set of units

			* metal:

			want kB in eV/K
			kB = 8.617e-5 eV/K

			* real:

			want kB in Kcal/mole/K
			kB = 0.001987191 ~= (8.617e-5 eV/K) (1 Kcal / 2.613e22 eV) (6.023 atoms/mole)

			(1) force to velocity conversion (in integrator)

			v = c (tf/m)
			c = force->ftm2v = converts ft/m to velocity

			* real:

			ft/m is Kcal/mole-Ang fmsec/(g/mole)
			v is Ang/fmsec

			ft/m (cgs) = ft/m (4.184E10 erg / Kcal) (sec / E15 fmsec) (E8 Ang / cm)
			v (cgs) = v (cm / E8 Ang) (E15 fmsec / sec)

			c = 4.184 E-4 = 1 / (48.88821 * 48.88821)

			* metal:

			ft/m is eV/Ang psec/(g/mole)
			v is Ang/psec

			ft/m (cgs) = ft/m (1.60219E-12 erg / eV) (sec / E12 psec)
			(6.022E23 atoms / mole) (E8 Ang / cm)
			v (cgs) = v (cm / E8 Ang) (E12 psec / sec)

			c = 9.648388E3 = 1 / 1.036443E-4

			(2) mv^2 to KE conversion (in thermo)

			boltz T = c (m v^2)
			c = force->mvv2e = converts mv^2 to energy

			* real:

			mv^2 is g/mole Ang^2/fmsec^2
			eng = Kcal/mole

			mv^2 (cgs) = mv^2 (mole / 6.022E23 atoms) (cm^2 / E16 Ang^2)
			(E30 fmsec^2 / sec^2)
			eng (cgs) = eng (4.184E10 erg / Kcal) (mole / 6.022E23 atoms)

			c = 2390.057 = 48.88821 * 48.888821

			* metal:

			mv^2 is g/mole Ang^2/psec^2
			eng = eV

			mv^2 (cgs) = mv^2 (mole / 6.022E23 atoms) (cm^2 / E16 Ang^2)
			(E24 psec^2 / sec^2)
			eng (cgs) = eng (1.60219E-12 erg / eV)

			c = 1.036443E-4

			(3) NkT/V to pressure conversion (in thermo and NPT integrator)

			P = c (NkT/V)
			c = force->nktv2p = converts NkT/V to pressure

			* real:

			NkT/V is Kcal/mole-K K/Ang^3
			P is atm

			NkT/V (cgs) = NkT/V (mole / 6.022E23 atoms) (4.184E10 erg / Kcal)
			(E24 Ang^3 / cm^3)
			P (cgs) = P (1.013E6 dyne/cm^2 / atm)

			c = 68586.95

			* metal:

			NkT/V is eV/K K/Ang^3
			P is bar

			NkT/V (cgs) = NkT/V (1.60219E-12 erg / eV) (E24 Ang^3 / cm^3)
			P (cgs) = P (E6 dyne/cm^2 / bar)

			c = 1.60219E6

			(4) q^2/r to energy conversion (in pair potential)

			E = c (q^2 / r)
			c = force->qqr2e = converts q^2/r to energy

			* real:

			q^2/r is e^2/Ang
			eng is Kcal/mole

			q^2/r (cgs) = q^2/r (E8 Ang / cm) (1.60219E-19 Coul / e)^2
			(statCoul / 3.3356E-10 Coul)^2
			E (cgs) = E (4.184E10 erg / Kcal) (mole / 6.022E23 atoms)

			(statCoul is cgs unit of charge)
			c = 332.0697

			* metal:

			q^2/r is e^2/Ang
			eng is eV

			q^2/r (cgs) = q^2/r (E8 Ang / cm) (1.60219E-19 Coul / e)^2
			(statCoul / 3.3356E-10 Coul)^2
			E (cgs) = E (1.60219E-12 erg / eV)

			(statCoul is cgs unit of charge)
			c = 14.40659

			(5) qE to force conversion (in fix efield)

			F = c (q E)
			c = force->qe2f = converts qE to force

			* real:

			qE is e Volt/Ang = e joule/Coulomb-Ang, since Volt = joule/Coulomb
			F is Kcal/mole-Ang

			qE (cgs) = qE (1.60219E-19 Coul / e) (statCoul / 3.3356E-10 Coul)
			(1.0E8 Ang / cm) (1.0E7 erg / joule) (3.3356E-10 Coul / statCoul)
			F (cgs) = F (4.184E10 erg / Kcal) (mole / 6.022E23 atoms) (1.0E8 Ang / cm)

			(statCoul is cgs unit of charge)
			c = 1.0E-3 * 1.60219E-19 * 6.022E23 / 4.184 = 23.0602

			* metal:

			qE is e Volt/Ang = e joule/Coulomb-Ang, since Volt = joule/Coulomb
			F is eV/Ang

			qE (cgs) = qE (1.60219E-19 Coul / e) (statCoul / 3.3356E-10 Coul)
			(1.0E8 Ang / cm) (1.0E7 erg / joule) (3.3356E-10 Coul / statCoul)
			F (cgs) = F (1.60219E-12 erg / eV) (1.0E8 Ang / cm)

			(statCoul is cgs unit of charge)
			c = 1.0E7 * 1.60219E-19 / 1.60219E-12 = 1.0

			--------------------------------------------
			1 Jun 2004
			Roy Pollock on non-neutral systems via PPPM

			>Hi Roy - quick question on charge-neutrality in PPPM, that someone is
			>asking about. Here is my understanding of what PPPM (or standard
			>Ewald) does - is this correct ?
			>
			>Thanks,
			>Steve
			>
			>-------------
			>
			>qsum = sum of user-specified charges
			>
			>* If the user specifies a charge-neutral system:
			>
			>the energy as computed in PPPM is:
			>
			> e_long = e_long - gewald*qsqsum/1.772453851D0 -
			> $ 0.5D0piqsumqsum/(gewaldgewaldxprdyprd*zprd_slab)
			>
			>and qsum = 0.0, so the last term is 0.0
			>
			>* If the user specifies a non-neutral system:
			>
			>then a compensating uniform background charge (opposite sign) is added
			>by LAMMPS, else the system energy blows up. Qsum is non-zero and the
			>last term in the above formula is the energy term due to the
			>compensating background charge.
			>
			>The forces on each atom are the same in both cases since a non-neutral
			>system has the same forces on each atom ...
			>
			Steve,

			Yes, Sounds correct to me. Since we're treating an infinite (periodic)
			system you need the charge
			neutralizing background.

			Treating either electrons (or ions) as a smeared uniform
			background is done in a lot of condensed matter or plasma physics
			models but not, I guess, in biological systems.

			Roy

			If you want a quick numerical check on PPPM for the uniform
			background case (OCP)
			a couple of static lattice Coulomb energies are:

			U/N=-1.41865 Z**2/r_nn/r_nn for a simple cubic lattice
			and
			U/N = -1.57583 Z**2/r_nn for a bcc lattice

			where r_nn is the nearest neighbor distance.

			If it gets that right for the uniform background case then its
			probably coded correctly.
			If not I'm going into the witness protection program. (I may do it anyway).

			Roy

			--------------------------------------------
			26 Apr 2004
			PPPM for dipoles with Roy

			I haven't permanently forgotten your inquiry about PPPM dipoles.
			I'm presently in over my head trying to finish off a paper on nuclear
			reaction rates in plasmas (more interesting than it sounds actually)
			but do want to eventually incorporate dipoles into the PPPM
			algorithm. The main theoretical bit is, probably, to derive the modified
			Green's function to compensate for using the non-Gaussian assignment
			function in assigning moments to grid points ?
			Hope this was a long term project.
			Roy

			told him I'd let him know how my 2-particle implementation went

			>So I presume I should have each dipole contribute 2 Q^2 to the summed
			>squared charge calculation, as it affects G_ewald ??

			Yes. The 2 charges, Q and -Q, are just 2 independent charges as far
			as the formula is concerned. If you analytically take the limit of
			the spacing going to zero you'll just rederive the Ewald (Korngold,
			Kornblatt, somebody?) formula for dipoles.

			You want the 2 charges representing a dipole to look like a pure
			dipole at the nearest neighbors. Just writing down the potential from
			the 2 charges and expanding in their separation, call it a, this
			requires a/Rnn<<1 (Rnn nearest neighbor distance). The charge Q
			doesn't affect the inequality based on the ratio ratio of dipole to
			higher multipoles in the expansion.

			--------------------------------------------
			5 Jan 2004
			John Carpenter ideas

			wants to create a LAMMPS users group

			It would be really good if we could get a LAMMPS User Group
			going. It could a rallying point for users to contribute
			software (i.e., useful diagnostic routines, other converter
			programs, dump file analysis routines, etc.). I would be
			willing to help in organizing such a Group and maybe even
			maintaining a FTP server for contributed software, etc.
			What do you think?

			Having a User Group doesn't necessarily put
			a demand on you, but a User Group that the
			code author doesn't participate in is probably
			stillborn. One might imagine a possible User
			Group meeting associated with a national
			scientific meeting as consisting of a short
			talk by you about the next version and maybe
			a few short talks by other users. It would
			mainly be a chance for users to meet one
			another and exchange ideas and possibly code.

			bug in 2001 with ndof in KE vs temp vs velocity creation

			In particular, in the test job class2/in.class2
			(which I don't have on my system - early version of
			LAMMPS 2001?) the input file specifies an initial
			temperature of 300, but the initial temperature
			(step 0) printed out comes up with 313 versus
			300 in the reference log file, but the kinetic
			energy for both agreed. We got them to agree by
			deleting the (-3) from the ndof calculation in
			thermo.f

			3M desire for COMPASS ff

			I've gone to 3M a couple of times and got them
			going. They had some trouble with msi2lmp and
			lmp2arc. They also had some misconceptions about
			a given atom interacting with images of itself
			outside the box, i.e., if you put a single water
			molecule inside a small box, you should get the energy
			of crystalline water. Incidentally, the minimizer
			core dumps on this. I haven't had a chance to
			track that down yet. Nevertheless, they are up and
			running now. They have plans to implement the
			COMPASS forcefield in LAMMPS. COMPASS is the Class II
			field developed at Cal Tech for the polymer sciences
			group at Accelys. 3M was dismayed that LAMMPS 2003 does
			not yet support Class II forcefields. I have a copy of
			the COMPASS paper and Sanat Mohanty at 3M and I will
			discuss the tasks that are needed to do that. There
			are two tricky parts: a msi2lmp equivalent and the
			H bond terms in the forcefield. Other than that the
			rest of the terms are all there.

			--------------------------------------------
			1 Jan 2004
			Crack simulations by Chris Kimmer at SNL/CA

			using MD - see Seminars

			--------------------------------------------
			1 Jan 2004
			PreWarp, MEAM comments from Greg Wagner

			I think the main tool you're refering to is Prewarp, which is a C++ code
			that's mainly the work of Phil Gullett. Phil's in the middle of moving
			to a faculty position at Mississippi State (aka "Sandia South"), but I'm
			cc'ing him in case he can help.

			I haven't spent much time looking at the Prewarp code, but since it's a
			C++ code it would probably be pretty easy to copy some classes from
			LAMMPS into Prewarp to handle output to a format that LAMMPS can read.
			Phil, is this reasonable?

			As for MEAM, I'm finding that I just don't use it much anymore, since I
			never did get a good set of Ni-H parameters working. It's too bad,
			since I did spend quite a bit of time on it. Have you had many requests
			for it in LAMMPS?

			I said:
			I recall you made some other additions to Warp, that might be good
			to put into the new LAMMPS - e.g. Ensight output format, etc. Don't know
			what they all were, however ...

			Anyone who's asked me for MEAM, I've referred to you. I thought the
			LANL folks used it quite a bit, but I haven't talked to Srini in many
			months. I guess we should wait until someone (SNL or outside) is
			motivated to add it themselves.

			--------------------------------------------
			1 Jan 2004
			Goddard work with LAMMPS and Dreiding ff

			Jang04.pdf paper

			from Paul:
			Steve, I was surprised to see that they are using LAMMPS with the
			Dreiding FF. I wonder if we could persuade them to share their "upgrade"
			of LAMMPS with us. It would be nice to include the Dreiding FF with
			LAMMPS. Did you know about this work?

			I assumed that the functional forms were different because they said
			they modified LAMMPS to handle their FF. Is the factor of 2 in the EvdW
			what they are talking about? That seems to be non-standard LJ potential.
			The table of coeffs for this particular system are included in the
			paper. It would be nice to have all of their parameters though. -Paul

			--------------------------------------------
			1 Jan 2004
			divergence with F90 GranFlow

			see emails from Gary - 13 Apr 04
			they have a few new features and attractive potential and analysis
			are working on lattice Boltzmann or pseudospectral in serial
			Jeremy Lechman (jblechm) is new guy

			So far as F90 vs C++, here is the status. The new version of LAMMPS
			is C++ and I've put many of the F90 GranFlow code features into it:
			granular potentials, boundary conditions, etc. However there are some
			features not yet working - like storing shear history in restart
			files, that I haven't had time to code up. Also, Gary indicated there
			are some new GranFlow features that I haven't seen - so it may take a
			bit of work to get the C++ version up to full functionality.

			I don't know how comfortable you are with C++ vs F90, so that's
			another issue.

			Having said that, I would still vote to add new functionality (like PS
			or lattice-Boltzmann) to the C++ version. I think it will be easier
			to add and more supported in the long run. But I'm open to discussing
			the trade-offs. One thing I would like to do is get you or Gray to
			use the granular features in the new LAMMPS so you can test/break
			things and give me feedback ...

			--------------------------------------------
			1 Jan 2004
			classes that create & manage their own atom-based arrays

			via callbacks (so are all nmax in length):

			fix_shake: shake_flag, shake_atom, shake_type, xshake
			fix_shear_history: npartner, partner, shear
			fix_rigid: displace
			fix_wiggle: original
			diag_diff: xoriginal

			within the class (so can be nlocal or nmax in length):

			verlet: f_pair[nmax]
			neighbor: bins[nmax], numneigh[nlocal], firstneigh[nlocal], xhold[nlocal]
			ewald & pppm: several arrays that are [nlocal] in length

			--------------------------------------------
			1 Jan 2004
			Overloaded classes that are derived from parents

			pair_lj_charmm_coul_charmm_implicit from pair_lj_charmm_coul_charmm
			pair_gran_no_history from pair_gran_history
			pair_gran_hertzian from pair_gran_history
			pair_lj_cut_coul_dh from pair_lj_cut_coul_cut

			--------------------------------------------
			1 Jan 2004
			Code locations where I have to know type of child class
			instead of just using generic parent

			in temper.cpp:
			set FixNVT or FixLangevin to get at reset_target()
			in temperature.cpp:
			set FixShake to get at dof()
			in pressure.cpp:
			set FixShake to get at virial array
			in fix_nve_dipole.cpp
			set PairDipole to get at sigma array
			in dihedral_charmm.cpp::init()
			set pair to one of 3 pair styles to grab lj14 coeffs
			in fix_wall_gran.cpp
			set PairGranHistory (and other 2) to get at xkk,xkkt
			in pair_gran_history.cpp
			set FixInsert to get at radius_hi
			in fix_insert.cpp
			access FixShearHistory fix to get at its npartner array
			in fix_gran_diag.cpp
			set PairGranHistory (and other 2) to get at xkk,xkkt,etc
			in neighbor.cpp
			set FixShearHistory to get at npartner,partner in pair_granular fns

			--------------------------------------------
			1 Jan 2004
			how bond types (angle, dihedral, improper) are set to 0 or negative
			for breaking, turning off, etc

			functions that alter bond type

			delete_bonds sets all of them negative
			if specify remove, it removes anything <= 0
			fix_shake sets bond,angle negative
			when class is deleted (unfix or exit), sets them back positive
			bond_quartic sets bonds to 0 when a bond breaks
			also skips computation of any bond that is <= 0
			neighbor tests if any bonds are <= 0 now or could be during run
			if so, calls bond_partial function
			bond_partial: if bond type <= 0, does not store bond in bond neigh list
			restart prints all bonds, but always writes out bondtype as >= 0
			so would have to reset the fix_shake or delete_bonds
			but deleted bond_quartic would still be in restart file
			dump bond prints out any bond with type != 0, including negative types

			--------------------------------------------
			1 Jan 2004
			hard-to-fix "bug" with SHAKE and NPT

			the SHAKE virial is not set until SHAKE is called the 1st time
			the setup in verlet, calls the setup for all fixes
			if NPT appears before SHAKE in the input script then
			NPT setup will call pressure which expects SHAKE virial to be set
			to keep from using non-init variable, I set shake virial = 0.0
			when shake fix is created
			if NPT appears after SHAKE in input script then
			all is correct
			SHAKE setup occurs 1st and virial is computed

			so only "bug" is that the initial pressure computed by NPT setup
			is slightly different depending on order of fixes in input script
			doesn't look different on timestep 0 thermo printout, but
			will probably affect dynamics slightly via NPT

			only way to correct this, would be to order the setup invocations
			depending on what fixes the user had defined (force NPT after SHAKE)
			don't really want to mess with dependencies like that

			--------------------------------------------
			1 Jan 2004
			use of temperature class

			used by 3 fixes (nvt, temp/rescale, npt)
			not by Langevin since it doesn't compute a temp
			fix_modify will change what temp ID is used in these commands
			warn if temp group != fix group
			when fix is created:
			if group is all, uses default temp
			else creates group with fictitious name (FixNVT + id + group-name)
			and temp style "full"
			if fictitious name already exists, just point to that one
			used by 2 velocity (create and scale)
			velocity args can change what temp ID is used
			if temp arg not specified, then velocity uses default temp
			but changes the group setting inside temp to those of the velocity group
			changes them back at end of velocity command
			used by thermo
			thermo_modify will change what temp ID is used

			when defined a temp class has both a style (full, partial, etc)
			and a group - is given a name by user
			default = style full and group all

			all temps are init once (from within force) when a simulation starts

			whenever a temp is init, it checks for a SHAKE fix and counts SHAKE dof
			only done in its group
			means that velocity must have SHAKE fix defined before it
			else will miss those dof when computing temp

			no easy way to delete a temperature:
			would have to repack the other ones
			this would change ptrs
			they are all deleted by force at LAMMPS exit

			--------------------------------------------
			1 Jan 2004
			General C++ questions (ask Roscoe):

			made atom,memory,etc as static members of globally-inherited System class
			this is so there will be only one copy and everyone can see it
			is this best way or only way to do it?
			makes lammps.cpp and system.cpp and system.h kind of ugly
			any other way to make top-level System visible in lammps.cpp

			when to use virtual classes versus just derived classes
			for nonbond, integrate, fix, dump, etc
			why not just have atom be a regular class with dummy defs for its
			functions (including some fns that are only in some child classes, like
			eam)

			doing neighbor styles with fn ptrs in Neighbor class
			any way to make a sub-class on the fly from within Neighbor,
			that inherits access to Neighbor data (e.g. pages)

			is it efficient/legal in read_data.cpp to say Atom &atom = *sys.atom;
			why can't put Atom &atom in read_data.h?
			can do it with pre-initializer in constructor is only way

			using a fn pointer (or class, e.g. pair->compute()) to replace
			a big if test in F90 - is that efficient in C++

			--------------------------------------------
			1 Jan 2004
			Version comparisons:

			any NPT & vol_rescale will not agree in pair energy, since f90 computes at end
			of timestep, C++ computes at middle (box size is different)
			does this mean long eng is different too?

			in.lj.nve matches for 1000, then diverges slowly
			in.lj.npt does not agree in energy (see above)
			in.lj.langevin - have to comment out fix_langevin::pre() call to force()
			to get agreement with F90 version
			in.bead - have to invoke fix_langevin::pre() call to force()
			to get agreement with F90 version
			in.lc.small.periodic - slight diff in pressure at step 100
			in.lc.small.nonperiodic - seem OK for 1000 steps
			in.stouch (with NVT) - slight diff in pressure at step 100
			in.stouch (with NPT) - seems to be drifting faster than would expect
			already different a bit by 10 timesteps
			in.polymer.chain.a - match to 1000 steps, different at 2000, diverges after
			in.polymer.chain.ab - starts to unmatch about 1000 steps
			in.lj.volume - does not agree in eng (see above) or when re-neigh
			in.lipid - OK for 1st few 100 steps - slower on Linux, OK on tflops

			--------------------------------------------
			1 Jan 2004
			Catching a new error:

			this will fail gracefully on huge malloc request

			#include "stdio.h"
			#include "stdlib.h"
			#include <new>
			//using namespace std;

			main (int argc, char **argv)
			{
			int n = atoi(argv[1]);
			printf("New on %d ints\n",n);

			int *a = NULL;
			try {
			a = new int[n];
			}
			catch (std::bad_alloc xa) {
			printf("ERROR: Failed malloc %p\n",a);
			}

			printf("Ptr to array = %p\n",a);
			// a[0] = 1;
			//a[n-1] = n;
			//printf("Array values %d %d\n",a[0],a[n-1]);
			delete [] a;
			}

			--------------------------------------------
			1 Jan 2004
			ideas for a presentation on new LAMMPS

			motivation:
			(1) working on LDRD to add dipoles
			what that affected w/in LAMMPS (integrator, new atom fields, comm, thermo)
			always worry when add something that
			a) something else will break
			b) code will slow down, become less parallel
			LAMMPS had global vars, reasonably modular, but big and unwieldy
			(2) maintain 3 codes: warp, LAMMPS, Gran - similar in basics, but
			all different atom styles/arays, force fields, BC, file formats
			LAMMPS 2001 = 34700 lines
			Warp = 11600 lines
			ParaDyn = 6200 lines
			GranFlow = 10100 lines
			(3) working on C++ in other projects, started to know enough to be dangerous

			C++ lets me combine features of all 3 codes in modular way:
			my C++ is C underneath, C++ at top to organize it
			basic code architecture:
			system class inherited by everyone
			contains neighbor, force, atom ptrs visible to all

			nice C++ (or C) features for a code like LAMMPS:
			might be ways to do some of these in F90, not saying is impossible
			some cases are possible, just aren't the way people write Fortran
			some cases is clunky in F90, or I just don't know how
			fn ptrs or class instances:
			big if test in integrator over pair potentials replaced by one statement
			pair->compute
			show code snippet from integrate routine in src vs new
			everything about a new option is in one place, one file
			in F90 all info was scattered thru code - input, start, fix
			dynamic memory - neighbor lists - realloc
			totally dynamic memory (should never run out in mid-run)
			multiple instances of things:
			several RNG in code, several temperatures, etc
			easy in C++, not so easy in F90
			a class can have one interface which rest of code knows
			underneath it can morph into lots of things
			this is inheritance, not possible in F90
			examples: atom, update, dump
			talk about what interface is to those
			an atom needs to know:
			how to read from file
			how to allocate/grow arrays for them
			how to pack/unpack for moving to new proc
			what to communicate
			a fix needs to know:
			what to do at various stages of timestep
			large # of MD options can be cast as a fix
			give examples: nve, npt, aveforce, gravity, wall, shake, bondswap
			integrator doesn't need to know how they work, when to call
			simplifies integrator greatly
			dump needs to know:
			how to write header
			how to package up
			how to write a line to a file (format, specific info)

			new features of LAMMPS you may like:
			a couple of these are "Grest mode" features
			heard Gary was submitting jobs from his hospital bed
			meant he was editing files, didn't want him to die in mid keystroke
			these features will let him re-submit jobs with less effort, editing
			read restart from file*
			cd, jump, clear
			groups: used for BC & fixes, output, init of velocity
			can have multiple NVT groups - e.g. a McDLT fix (hot/cold on 2 sides)
			can make several output dumps of different atoms and different
			quantities to different files
			variable (2 kinds), next, proc clusters
			show simple examples
			will let you run 8 jobs of 16-each procs on 128 procs
			and run 100 jobs on those 8 clusters, one after the other
			will let you grab 128 procs for one week from now, fill in
			low-level script later
			tempering (with M Sears)
			tell how it works
			show plots of overlapping PE
			how it works on a peptide or protein complex, Garcia-like results
			have ability to take PDB file and solvate it, convert into LAMMPS input
			xmovie - available if you want it
			little animator in Qt - available if you want it
			not competitor to Rasmol, Raster3d, or VMD
			post-processing converts also available for those
			rigid fixes

			advantages of LAMMPS in this mode:
			easy for me or user to add to without impacting anything else
			don't have to know C++ to do it, just use another as template, write C
			things that can easily be added:
			force field
			atom style
			new integrator
			new thermo (temp, press)
			new dump
			new BC (fix, ensemble)
			new diagnostic (is it a fix?)
			one code base, code can do all that other 3 codes did
			illustrate w/ simple movie of 3-4 kinds of calculations
			size of new code versus sum of old 3
			could add hybrid force field (EAM + CHARMM) - tell how

			is it as fast as before:
			show some timing curves
			serial and parallel
			low-level routine is C instead of F90
			kept all LAMMPS speed-effects
			neighbor lists
			ghost atoms
			6-way communication (cuts corners, still true even for long cutoffs)

			what it can't do:
			other parallel styles of comm (force, atom) (ParaDyn)
			load-balancing
			reactive FF
			reacXX is a very complex "pair" potential
			charge equil
			talk about how similar/different to GRASP
			not doing multiple force fields & neighbor lists on sub-groups of atoms

			what I miss from F90:
			multi-dim arrays - e.g. allocate(f(3,n),bonds(2,3,n))
			cs(-4:4,3,1000)

			-----------------------------------------
			1 Jan 2004
			Class overview

			0. LAMMPS

			1. system - all next-level classes inherit from it
			static ptr to single instance of these:
			atom, modify, force, update, neighbor, comm, domain,
			group, output, timer, error, memory
			these ptrs are initialized within lammps.cpp

			1a. atom - contains all atom-based data

			1b. update - integrators/minimizers
			integrate - one of these
			minimize - one of these

			1b1. integrate:
			virtual class, derived ones are Verlet, respa

			1b2. minimize:
			virtual class, derived ones are CG, HTFN

			1c. force - all force and energy routines
			nonbond - one of these
			bond - one of these

			1c1. nonbond - nonbond force routine
			virtual class, derived ones are LJ, soft, etc

			1c2. bond - bond force routine
			virtual class, derived ones are FENE, harmonic, etc

			1d. modify - ensemble, fix routines
			fix - array of these

			1d1. fix - any routine that changes v,f,x
			virtual class, derived ones are NVE, NPT, umbrella, etc

			1e. neighbor - neighbor list routines

			1f. comm - communication routines

			1g. domain - physical simlation box

			1h. group - defined groups of atoms

			1i. output - all output to screen/files
			dump - array of these
			thermo - one of these

			1i1. dump - output of vel, atom, force
			virtual class, derived ones are atom, velocity, force, etc
			generic write function in virtual class
			write_header, pack, write_data are specific to derived classes

			1i2. thermo - thermodynamic computations
			temperature - one of these
			pressure - one of these

			1i2a. temperature - compute temp, need to be virtual?
			1i2b. pressure - compute pressure, need to be virtual?

			1j. timer - CPU timings

			1k. error - error handlers

			1l. memory - memory allocation routines

			2. input - grab and parse lines from input script
			created/deleted in lammps.cpp, used nowhere else
			returns next command to execute

			3. read_data - reads input data file
			created/deleted in lammps.cpp

			4. velocity - sets velocity of some atoms
			created/deleted in lammps.cpp

			5. random - RNG
			created/deleted within Velocity class

			6. special - create molecular topology
			created/deleted within ReadData class

			-----------------------------------------
			1 Jan 2004
			General documentation

			How to run without (potentially) gigantic map() array in atom.cpp
			only OK to do for LJ or EAM problems that are non-molecular
			and don't use map() in special or SHAKE or ??

			in atom.cpp, are 3 routines: map_init(), map_clear(), map_set()
			comment out the code in each of these 3 routines
			then LAMMPS will not allocate this array

			My question on bond frequencies:

			> I'm looking for bond frequency/periods - I found a table that lists
			> their absorption spectra in units of inverse-cm. E.g. C-H = 2900
			> 1/cm, O-H = 3000 1/cm.
			>
			> Do you know the conversion to time units? E.g. I want 10 fmsec for
			> the characteristic bond period.

			Answer, from Paul:

			Just multiply cm^(-1) by c to get per-sec = frequency

			e.g. 1/(310^10 cm/s 3000 cm^-1) = 11.11 fs

			Also see:
			http://www.chem.umd.edu/courses/manoso/Lectures/IR.pdf

			C-O = 2143 cm-1
			NO = 1875 cm-1
			CC = 1200 cm-1

			------------
			Documentation for screen and logfile settings

			depends on -screen and -log command-line args, 1 vs multi world,
			which proc I am in universe and local world

			uscreen = universe::uscreen
			screen = system::screen
			ulogfile = universe::ulogfile
			logfile = system::logfile
			N = which world (0 to nworld-1)

			universe has 1 world:

			universe proc 0
			uscreen = stdout
			if -screen, uscreen = file (NULL if none)
			ulogfile = log.lammps
			if -log, ulogfile = file (NULL if none)
			universe proc 1-N
			uscreen = stdout
			if -screen, uscreen = NULL
			ulogfile = NULL
			world proc 0
			screen = stdout
			if -screen, screen = file (NULL if none)
			logfile = log.lammps
			if -log, logfile = file (NULL if none)
			world proc 1-N
			screen = stdout
			if -screen, screen = NULL
			logfile = NULL

			multi-world universe:

			universe proc 0
			uscreen = stdout
			if -screen, uscreen = file (NULL if none)
			ulogfile = log.lammps
			if -log, ulogfile = file (NULL if none)
			universe proc 1-N
			uscreen = stdout
			if -screen, uscreen = NULL
			ulogfile = NULL
			world proc 0
			screen = screen.N
			if -screen, screen = file.N (NULL if none)
			logfile = log.lammps.N
			if -log, logfile = file.N (NULL if none)
			world proc 1-N
			screen = NULL
			logfile = NULL

			the "log" command simply closes logfile on each world 0
			reopens the system::logfile with the specific name only on world proc 0
			does not change universe::ulogfile

			--------------------------
			1 Jan 2004
			NAMD benchmark

			3 bonds per water makes SHAKE not work with "m" switch
			NVE vs NVT
			reneigh criterion
			64^3 mesh is too hi precision
			run on Lemeiux, is 2002 Lem same as now?
			set temp to not include SHAKE so get accurate temp at initial step

			-----------------------
			1 Jan 2004
			Aidan ideas

			multiple neighbor lists (for diff force fields, cutoffs)
			only applied to some atoms
			could I do that for REPSA short range (bond only)
			also need diff comm for diff neighbor distances
			have multiple force fields (a, b, hybrid ab)
			also only applied to some atoms

lammps/src/Notes284c1f2d04aclammps-icms

NotesView Options

lammps/src/Notes
284c1f2d04aclammps-icms

Notes
View Options