Quasicontinuum Code Revision History
Version 1.4 (November 2011)
A number of bugs were fixed:
- Fixed a bug that could cause some crystals to be built incorrectly
when the orientations of the global axes are along high-index crystal
- Corrected a bug that could cause the connectivity of repatoms to
elements to be incorrectly computed for some rather special mesh
- Corrected a small bug in the calculation of the load-displacement
curve for the nano-indentation example.
- Corrected parsing of input lines in qcmain.f. Blank lines are
no longer placed in qc.cmd
- Corrected parsing of input lines in mod_head.f. Old code could
crash if blank lines appeared in qc.cmd
- Allocation of grad and bsave corrected in mod_solve.f. Old code
would crash if "solv,nr" was called after a "solv,cg".
- A bug in mod_repatom was corrected that results in a segmentation fault
whenever the size of arrays used to find neighboring elements in
BuildSeenElementsList has to be extended.
- A bug in mod_repatom was corrected that occasionally led to missed
atomic neighbors when there are multiple grains.
- Bug in mod_plotting fixed. The old code did not correctly generate
strain patch plots (no strains appeared in the file).
- Two new features were added:
- Added a macro "mark" that sets or returns to a marked position where the
displacement field is stored.
- Added a convergence test on the number of nonlocal nodes
('nnln'). This item can be used together with the previous item to back-up in the event of a nucleation event and to restart from the displacement field prior to the event.
Several optimizations were made:
The routine to determine in a point lies inside a given polygon (PointInPoly) was optimized.
The routine that builds the connectivity array storing all elements touching a given repatom (BuildIrel) was optimized.
The adaption routine were improved so that the refinement of elements on which boundary conditions are applied is more robust (improved routine "fixcon").
Version 1.3 (May 2007)
The QC Code was largely rewritten and restructured to conform to the
Fortran 90 standard. As a result the code is very portable. It has been
tested with ten different Fortran compilers. In all cases it compiles and
runs without modifications. Some of the main changes introduced are:
- All functions and subroutines were placed within modules. This ensures
explicit interfaces for called subroutines and functions.. As a result,
a variety of potentially non-portable Fortran 77 "tricks" no longer work.
These have been removed.
- "double precision" variables were replaced with "real(kind=dp)", where
dp is defined in mod_global as "selected_real_kind(p=15,r=307)". This
ensures that the same precision is used on any machine. The change
to "dp" means generic routines must be used such as sin(x), sqrt(x)
and not the specific double precision routines dsin(x), dsqrt(x), etc.
Also, type change to "dp" precision is done with real(n,dp), where n is
an integer, and not dble(n).
- Obsolete features in Fortran 90 were removed. For example,
"character*n" were changed to "character(len=n)". Arithmetic
if statements removed, etc.
- Passing of array sections to subroutines were replaced by passing a
temporary array containing the relevant section. Technically, this
is not necessary since Fortran 90 allows passing array sections.
However, in practice some compilers failed to do do this correctly
especially, when the array was part of a derived type.
- Dynamic memory allocation of variables was done automatically
whenever possible, otherwise the "allocatable" clause was used.
The use of pointers was limited to cases where this was unavoidable.
Greatly improved streamlined make architecture designed to allow easy
porting from one platform to another. The new approach is very easy
to use and robust. It has been tested on a variety of Unix and Linux
systems as well as OS X 10.4.9.
The code was optimized in a number of ways. A speed-up of about 30-40%
was observed for typical applications when switching from version 1.2 to
The memory requirement of the code were greatly reduced. This was done
by doing the following:
- A number of memory leaks were detected and fixed.
- The irel() array was reduced from storing all neighbors of
a repatom within the cutoff radius to only its nearest neighbors.
- The arrays xl() and ul() were removed.
Restart file size was greatly reduced due to the changes in item 4. Note
that version 1.3 will still read in restart files from older versions
of the code. This is done automatically without user intervention.
Numerous bugs fixed. Most are minor and will not be mentioned. The more
major ones are listed below. The subroutine containing the bug appears
- (StatusCalc) The calculation of deformation gradients and
eigenvalues could be corrupted after adaption. As a result
repatom status (local or nonlocal) could be incorrectly computed.
This would be corrected in subsequent steps but could explain some
"jerky" behavior that was sometimes observed in simulations. This was
fixed. Note: The corrected code gives results that can be different
from earlier versions because the set of nonlocal atoms can be different.
- (pmacr) When calling restart b() was stored in blstat()
for 'read' but not for 'txrd'. This was fixed.
- (delaunay) Bug related to convex hull triangulation for the special
case where ncb=0 and nce/=0 fixed.
- (GetCellData) For "bad" crystallographic orientations with very large
repeat distance the algorithm to find the repeating cell could fail.
This was fixed.
A number of changes were introduced to the structure of the required
user routines that interface with QC. The changes make the routines
easier to work with and safer (i.e. less likely to "break" QC due
to any errors in them). The main changes are:
A file user_template.f is included in the QC directory as an aid to
porting user routines from older versions to version 1.3. The file
contains a skeleton of the file along with instructions on where and
how to place the older user routines.
- All user routines are now in a single file called "user_APP.f" where
APP is a particular application (e.g. user_gb.f). This file contains:
user_mesh -- mesh generator
user_bcon -- boundary conditions
user_pdel -- force versus displacement
user_potential -- user-specified external potential
user_plot -- user-specified plot option
The last two are new in this release. See the reference manual for
- Direct access to many module variables is no longer possible to
prevent a user inadvertently "breaking" the code. Instead routines
were added to allow a user to query a module for the value of a
particular variable. These routines include:
NumGrains - return number of grains
GetGrainCellSize - Get the repeating cell dimensions for a given grain
GetGrainNumCell - Get the number of atoms in the repeating cell for a given grain
GetGrainCellAtom - Get coordinates of a given atom in the repeating cell for a given grain
GetGrainRefStiff - reference stiffness matrix for a given grain
GetGrainBvec - Get Bravais vector matrix for a given grain
GetGrainBinv - Get inverse Bravais vector matrix for a given grain
GetGrainRefatom - Get reference atom for a given grain
GetGrainSpecies - Get species for a given grain
GetGrainNumVrts - Get number of vertices for polygon of a given grain
GetGrainVertex - Get given vertex of the polygon of a given grain
GetGrainXlatvect - Get crystallographic axes for a given grain
GetLoadTime - Get current loading time
GetLoadTimeStep - Get current loading timestep
GetLoadPropFact - Get current loading proportional factor
GetLoadOldPropFact - Get previous loading proportional factor
- "only" clause added to all use statements.
- Calling header to subroutine qsortr() (called from user_mesh)
was changed to:
The variables 'ndim' and 'ind' that appeared between 'xkey' and 'sign'
in the old header were dropped.
The version 1.3 release is focused more on a restructuring of the code
as a basis for future development. However, there are a number of
minor extensions in this version:
Version 1.3 is downward compatible being able to read input files
used with earlier versions of the code.
- Two user routines user_potential and user_plot were added, allowing
the user to define a user-specified external potential and plot option.
See the reference manual for more information.
- The zone macro (formerly used to defined no-adaption zones) was
generalized to describe "special-attribute" zones. This currently
includes nonlocal zones and no-adaption zones. See the reference manual
- More input verification was added to the code to catch user input
The QC Reference Manual and QC Tutorial Guide were significantly
expanded. The Reference Manual, in particular, has an entirely new
chapter on the structure of the QC code and on writing the user
The utility program "dynpots.f" was updated from f77 to conform to
the f90 standard.
Version 1.2 (February 2005)
- Support for multiple unconnected meshed domains (This was included in
version 1.1 - but the new version is much more robust). A new example
"Friction-example" is included to demonstrate this capability. (RM)
- mod_crystal: Small bug in CryCount that would not have affected
performance was corrected (ET).
- New flag DummiesAllowed added. When set to .false. no dummy atoms are
allowed, meaning that all neighbors of nonlocal atoms must be repatoms
themselves. If the flag is set to false and dummies are detected in
a status calculation, adaption is triggered until there are no more
dummies. If using the NR solver no changes are necessary to the input
file besides setting the appropriate flag. For the CG solver the
solver must be placed inside an additional loop with a status,update
The reason is that neighbor lists are updated inside the CG solver.
This can result in dummy atoms being generated. It this happens and
dummies are not allowed, the CG solver will exit with a warning.
The loop above handles this through the status,update command
which will trigger a full status calculation and refine the mesh
to remove the dummies. (ET)
- Before writing out a restart file, a backup of the last written
restart file is saved (with "_LAST" added to its name). For example,
will copy the file qc_simulation.res (if it exists) to
qc_simulation_LAST.res, before writing out the new restart file.
The reason for this is that if something interesting happens in
a simulation, it is often helpful to have the restart file of the
step before the event occurred, in order to restart from there
with some changes to the input file, e.g. more adaption, etc. (ET)
- A new file called nonstandard.f has been added to the QC code.
This file contains any nonstandard Fortran extensions used by the
QC code to facilitate portability. Compilers which support
these extensions differently from most or that need to use specific
modules in order to access these extensions have their own nonstandard.f
file (nonstandard_compiler.f) along with a corresponding make file.
If you are porting QC to a new platform, modify nonstandard.f to support
the new compiler. (ET)
Currently nonstandard.f contains the following routines:
- flush command.
Flush is a non-standard Fortran command that forces the output buffers
to be written to disk. Without it the output files are only written when
the buffers become full (the size of the buffers is system dependent).
Because the command is non-standard, some compilers implement it
differently than others. So far we have encountered
trouble only with the SGI MIPSPro Fortran 90 compiler (versions 7.4
and higher) where the number of arguments to flush was changed.
As a result the previous version of the QC code crashed with a bus error.
(Thanks to Rafael Estevez at GEMPPM INSA Lyon for helping us figure this
- system command
System is a non-standard Fortran command that lets the program issue
an operating system command while running. The system command is
used to make the backup copy of the restart file described above.
- Various code optimizations were performed resulting in about a
25% speed-up (depending on the platform) (ET). These optimizations include:
- The function dot11 which was called by the active column solver
was removed and replaced with faster inline code.
- The function trarea which computes the area of a triangle was
eliminated and replaced with faster inline code.
- The function intri which determines whether a point is inside
a triangle was optimized. The part of the routine that
computed area coordinates and determined whether the point
was on the edges of the triangle was removed because this
information was never used.
- The routine fe_locate which finds out which element a given
point lies in was optimized.
- Some service routines in qclib (Dist2, distance, etc.) were
removed and replaced with faster inline code.
- The initial output generated by QC has been extended to include
the stress tensor and elastic constants matrix (in Voigt notation)
in the reference configuration. The stress tensor should be zero.
The degree it differs from zero is a measure of the accuracy of the
lattice parameter given to the program. Both the stress and elastic
constant components are given in relation to the global c.s. (ET)
To prevent the code from crashing in the event that moduli calculations
are not implemented in the potential being used, an internal error
code was added to the routines Calc_Local and NonLocal and to routines
calling them. This is transparent to the user and only affects developers.
- Improved reliability in neighbor finding for nonlocal repatoms.
In circumstances of very large relative deformation between two
neighboring grains or for grains with a highly non-convex shape
(for example with a long thin notch) the code could sometimes
miss neighbors to some repatoms, leading to an error in the
energy and force calculations. This has been fixed and the new
version is now extremely robust at the expense of some slowing
down of the code. (RM)
- If a model contained grains with different periodicities in the
out of plane direction, the code would not compute the neighbor
lists of nonlocal atoms correctly. The new version can now
handle this situation, although it is still constrained by the
2.5-D nature of the code and results must be interpreted with
- Multiple plots of the same data, for example the displacements
after each load step, can now be appended to multiple zones of a
single file, rather than multiple individual files. This is
achieved with the "append" flag in the "plot" command. (RM)
- The z-component of all repatom positions are now mapped onto Bravais
lattice sites within the first periodic copy of the crystal
structure. This doesn't change the results, but makes
visualization of out-of-plane deformations easier. (RM)
- The "nonlocal fully-refined" feature of the
adaption routine is now more robust in that it is less likely to
find circumstances where it cannot achieve the full refinement of
nonlocal repatoms. In the event that it fails, it terminates
with an error message. (RM)
- The input file is now much more flexible in terms of accepting
spaces at the start of commands and comment lines (denoted with a
"%") are now accepted everywhere in the file. (RM)
- A new option, 'type', has been added to the 'plot' command that
plots the positions and atomic number of each repatom. (RM)
Version 1.1 (May 2004)
- qclib.f: subroutine fe_locate was rewritten to allow for
multiple separate meshed regions. Previous versions were limited to a
single connected domain. (RM)
- mod_pload.f: The proportional load table was modified to handle a
step function loading without crashing (two prop values that have the same
time value). (RM)
- adaptgr2.f: Function GradientNorm() was modified so that an element
with a nonlocal node is automatically targetted for refinement even if
the centroid of the element lies inside a no-adaption zone. This is
necessary to prevent nonlocal nodes at the edge of no-adaption zones
representing more than themselves and generating an internal QC error.
(corrected by ET following a bug report by Markus Beuhler)
- mod_grain.f: (1) Target attribute added to grains so that pointers
to grain data can be defined later in subroutine nearestbravais in
polygrain.f. (2) Components of reciprocal lattice in global c.s.
computed and stored for each grain to save computation time later
in subroutine nearestbravais in polygrain.f.
(changes introduced by Jon Wilkening, Courant Institute of Mathematical
- mod_repatom.f: Data structures used by subroutine FindNeighbors to store
Bravais indices visited in neighbor search and dummy atoms were changed
from linked lists to pointers into vectors. The resulting code is
significantly faster. An additional speed-up was obtained by improving
the initial guess for the element breadth search used to locate the
element containing a specific atom. These changes will speed up status
calculations, neighbor list recalculations triggered in the solver and
(changes introduced by Jon Wilkening, Courant Institute of Mathematical
(1) The new initial guess used by the element search routine (fe_locate)
can result in slightly different results from those obtained using
previous versions of the code. The reason for this is that nodes lying
on the interface between two elements will sometimes be assigned to a
different element that in previous versions. As a result, due to
numerical precision issues, the energy and forces can be slightly
different. Because the solution procedure in the code is iterative,
even very small changes in forces can result in more pronounced changes
after many iterations. Qualitatively, the results are expected to be
the same, but small quantitative differences may be expected.
(2) A separate mod_repatom.f file for the Portland Group Compiler (PGF)
(mod_repatom.f.pgf) is no longer necessary. The new mod_repatom.f
file should work for PGF as well.
- "character*xx" statements were replaced with "character(len=xx)"
throughout the code, to prevent warnings in some fussy compilers. (ET)
- mod_poten: Bug corrected at end of shell loop in subroutine Calc_Local.
In the old version an incorrect error message was given in the rare
case where rinf=rcryst and the simulation was terminated. Also the
if statement checking whether to continue the loop did not work
correctly at the limit where all shells were visited. Both issues
were corrected. (ET)
- mod_solve: a "call output" was added in subroutine cgstep to inform
the user when the neighbor lists were updated during a CG step.
The updates were carried out in the previous version, but no message
was printed in the output stream. (ET)
- mod_poten: a bug which caused the code to crash when multiple
potential files were read in was fixed. The problem had to do with
an incorrect print statement. (ET)
- mod_poten: error handling was added so that if an incorrect
potential format was specified the code would exit gracefully with
a meaningful message. (ET)
- refine.f: this file containing a duplicate definition of function
refine was removed. (ET)
- plotting.f: A time and property stamp were addded to the title of Tecplot
plot files generated by the plot macro. This allows the user to know
which load step he or she is looking at from within Tecplot. (ET)
- Makefiles added for the Intel Fortran Compiler version 8.0. (ET)
Version 1.01 (August 2003)
- mod_repatomb.f: a minor bug was corrected (incorrect loop upper bound
on line 379).
- mod_adaptgr2.f: a minor bug was corrected (inconsistency in the
number of arguments in subroutine constraints on lines 248,252).
- qcmain.f: a minor bug was corrected (an extraneous comma in
the QC banner print statement was removed on line 157).
- A bug in the input file for the GB example (gb_shear.in) was corrected.
The crystallite radii for the grains was changed from 3.0 to 20.0
(lines 42 and 58).
- Make files for the Portland Group F90 compiler were added.
Note: The Portland Group F90 compiler generates error messages
when compiling mod_repatomb.f. This is because it does
not recognize the NULL() statement used to nullify pointers.
An alternative routine mod_repatomb.f.pgf has been included
in the distribution which solves this problem. If you are using
the Portland Group compiler you should replace mod_repatomb.f
with this file.
Version 1.0 (February 2003)