From a583fc1088a86867da7432fcb12c42f6a005639c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 8 Feb 2007 16:47:32 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@277 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/Notes | 2954 ----------------------------------------------------- 1 file changed, 2954 deletions(-) delete mode 100644 src/Notes diff --git a/src/Notes b/src/Notes deleted file mode 100644 index 71431f5ab1..0000000000 --- a/src/Notes +++ /dev/null @@ -1,2954 +0,0 @@ -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) 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" -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 + VMASS*VLOGV*VLOGV - (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.2*0.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 -#include -#include - -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 , 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.5D0*pi*qsum*qsum/(gewald*gewald*xprd*yprd*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 -//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/(3*10^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