ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into lammps-icms

Browse Source 
Resolved Conflicts:
	src/REPLICA/verlet_split.cpp
	src/USER-OMP/fix_omp.cpp
This commit is contained in:
Axel Kohlmeyer
2012-10-26 09:50:01 -04:00
parent 3930d66dfb 8663c92293
commit 5a75a6c827
25 changed files with 1511 additions and 1821 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
examples/COUPLE/README
View File

@ -32,7 +32,7 @@ lammps_spparks grain-growth Monte Carlo with strain via MD,
coupling to SPPARKS kinetic MC code
library collection of useful inter-code communication routines
simple simple example of driver code calling LAMMPS as library
fortran a wrapper on the LAMMPS library API that
fortran a simple wrapper on the LAMMPS library API that
can be called from Fortran
fortran2 a more sophisticated wrapper on the LAMMPS library API that
can be called from Fortran

2110
examples/COUPLE/fortran2/LAMMPS.F90
View File

File diff suppressed because it is too large Load Diff

486
examples/COUPLE/fortran2/README
View File

@ -1,221 +1,265 @@
LAMMPS.F90 defines a Fortran 2003 module, LAMMPS, which wraps all functions in
src/library.h so they can be used directly from Fortran-encoded programs.
All functions in src/library.h that use and/or return C-style pointers have
Fortran wrapper functions that use Fortran-style arrays, pointers, and
strings; all C-style memory management is handled internally with no user
intervention.
This interface was created by Karl Hammond who you can contact with
questions:
Karl D. Hammond
University of Tennessee, Knoxville
karlh at ugcs.caltech.edu
karlh at utk.edu
-------------------------------------
--COMPILATION--
First, be advised that mixed-language programming is not trivial. It requires
you to link in the required libraries of all languages you use (in this case,
those for Fortran, C, and C++), as well as any other libraries required.
You are also advised to read the --USE-- section below before trying to
compile.
The following steps will work to compile this module (replace ${LAMMPS_SRC}
with the path to your LAMMPS source directory):
(1) Compile LAMMPS as a static library. Call the resulting file ${LAMMPS_LIB},
which will have an actual name lake liblmp_openmpi.a. If compiling
using the MPI stubs in ${LAMMPS_SRC}/STUBS, you will need to know where
libmpi.a is as well (I'll call it ${MPI_STUBS} hereafter)
(2) Copy said library to your Fortran program's source directory or include
its location in a -L${LAMMPS_SRC} flag to your compiler.
(3) Compile (but don't link!) LAMMPS.F90. Example:
mpif90 -c LAMMPS.f90
OR
gfortran -c LAMMPS.F90
Copy the LAMMPS.o and lammps.mod (or whatever your compiler calls module
files) to your Fortran program's source directory.
NOTE: you may get a warning such as,
subroutine lammps_open_wrapper (argc, argv, communicator, ptr) &
Variable 'communicator' at (1) is a parameter to the BIND(C)
procedure 'lammps_open_wrapper' but may not be C interoperable
This is normal (see --IMPLEMENTATION NOTES--).
(4) Compile (but don't link) LAMMPS-wrapper.cpp. You will need its header
file as well. You will have to provide the locations of LAMMPS's
header files. For example,
mpicxx -c -I${LAMMPS_SRC} LAMMPS-wrapper.cpp
OR
g++ -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp
OR
icpc -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp
Copy the resulting object file LAMMPS-wrapper.o to your Fortran program's
source directory.
(4b) OPTIONAL: Make a library so you can carry around two files instead of
three. Example:
ar rs liblammps_fortran.a LAMMPS.o LAMMPS-wrapper.o
This will create the file liblammps_fortran.a that you can use in place
of "LAMMPS.o LAMMPS-wrapper.o" in part (6). Note that you will still
need to have the .mod file from part (3).
It is also possible to add LAMMPS.o and LAMMPS-wrapper.o into the
LAMMPS library (e.g., liblmp_openmpi.a) instead of creating a separate
library, like so:
ar rs ${LAMMPS_LIB} LAMMPS.o LAMMPS-wrapper.o
In this case, you can now use the Fortran wrapper functions as if they
were part of the usual LAMMPS library interface (if you have the module
file visible to the compiler, that is).
(5) Compile your Fortran program. Example:
mpif90 -c myfreeformatfile.f90
mpif90 -c myfixedformatfile.f
OR
gfortran -c myfreeformatfile.f90
gfortran -c myfixedformatfile.f
The object files generated by these steps are collectively referred to
as ${my_object_files} in the next step(s).
IMPORTANT: If the Fortran module from part (3) is not in the current
directory or in one searched by the compiler for module files, you will
need to include that location via the -I flag to the compiler.
(6) Link everything together, including any libraries needed by LAMMPS (such
as the C++ standard library, the C math library, the JPEG library, fftw,
etc.) For example,
mpif90 LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
${LAMMPS_LIB} -lstdc++ -lm
OR
gfortran LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
${LAMMPS_LIB} ${MPI_STUBS} -lstdc++ -lm
OR
ifort LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
${LAMMPS_LIB} ${MPI_STUBS} -cxxlib -limf -lm
Any other required libraries (e.g. -ljpeg, -lfftw) should be added to
the end of this line.
You should now have a working executable.
Steps 3 and 4 above are accomplished, possibly after some modifications to
the makefile, by make using the attached makefile.
-------------------------------------
--USAGE--
To use this API, your program unit (PROGRAM/SUBROUTINE/FUNCTION/MODULE/etc.)
should look something like this:
program call_lammps
use LAMMPS
! Other modules, etc.
implicit none
type (lammps_instance) :: lmp ! This is a pointer to your LAMMPS instance
double precision :: fix
double precision, dimension(:), allocatable :: fix2
! Rest of declarations
call lammps_open_no_mpi ('lmp -in /dev/null -screen out.lammps',lmp)
! Set up rest of program here
call lammps_file (lmp, 'in.example')
call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1)
call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1)
call lammps_close (lmp)
end program call_lammps
Important notes:
* All arguments which are char* variables in library.cpp are character (len=*)
variables here. For example,
call lammps_command (lmp, 'units metal')
will work as expected.
* The public functions (the only ones you can use) have interfaces as
described in the comments at the top of LAMMPS.F90. They are not always
the same as those in library.h, since C strings are replaced by Fortran
strings and the like.
* The module attempts to check whether you have done something stupid (such
as assign a 2D array to a scalar), but it's not perfect. For example, the
command
call lammps_extract_global (nlocal, ptr, 'nlocal')
will give nlocal correctly if nlocal is of type INTEGER, but it will give
the wrong answer if nlocal is of type REAL or DOUBLE PRECISION. This is a
feature of the (void*) type cast in library.cpp. There is no way I can
check this for you!
* You are allowed to use REAL or DOUBLE PRECISION floating-point numbers.
All LAMMPS data (which are of type REAL(C_double)) are rounded off if
placed in single precision variables. It is tacitly assumed that NO C++
variables are of type float; everything is int or double (since this is
all library.cpp currently handles).
* An example of a complete program is offered at the end of this file.
-------------------------------------
--TROUBLESHOOTING--
Compile-time errors probably indicate that your compiler is not new enough to
support Fortran 2003 features. For example, GCC 4.1.2 will not compile this
module, but GCC 4.4.0 will.
If your compiler balks at 'use, intrinsic :: ISO_C_binding,' try removing the
intrinsic part so it looks like an ordinary module. However, it is likely
that such a compiler will also have problems with everything else in the
file as well.
If you get a segfault as soon as the lammps_open call is made, check that you
compiled your program AND LAMMPS-header.cpp using the same MPI headers. Using
the stubs for one and the actual MPI library for the other will cause major
problems.
If you find run-time errors, please pass them along via the LAMMPS Users
mailing list. Please provide a minimal working example along with the names
and versions of the compilers you are using. Please make sure the error is
repeatable and is in MY code, not yours (generating a minimal working example
will usually ensure this anyway).
-------------------------------------
--IMPLEMENTATION NOTES--
The Fortran procedures have the same names as the C procedures, and
their purpose is the same, but they may take different arguments. Here are
some of the important differences:
* lammps_open and lammps_open_no_mpi take a string instead of argc and
argv. This is necessary because C and C++ have a very different way
of treating strings than Fortran.
* All C++ functions that accept char* pointers now accept Fortran-style
strings within this interface instead.
* All of the lammps_extract_[something] functions, which return void*
C-style pointers, have been replaced by generic subroutines that return
Fortran variables (which may be arrays). The first argument houses the
variable to be returned; all other arguments are identical except as
stipulated above. Note that it is not possible to declare generic
functions that are selected based solely on the type/kind/rank (TKR)
signature of the return value, only based on the TKR of the arguments.
* The SHAPE of the first argument to lammps_extract_[something] is checked
against the "shape" of the C array (e.g., double vs. double* vs. double**).
Calling a subroutine with arguments of inappropriate rank will result in an
error at run time.
* All arrays passed to subroutines must be ALLOCATABLE and are REALLOCATED
to fit the shape of the array LAMMPS will be returning.
* The indices i and j in lammps_extract_fix are used the same way they
are in f_ID[i][j] references in LAMMPS (i.e., starting from 1). This is
different than the way library.cpp uses these numbers, but is more
consistent with the way arrays are accessed in LAMMPS and in Fortran.
* The char* pointer normally returned by lammps_command is thrown away
in this version; note also that lammps_command is now a subroutine
instead of a function.
* The pointer to LAMMPS itself is of type(lammps_instance), which is itself
a synonym for type(C_ptr), part of ISO_C_BINDING. Type (C_ptr) is
C's void* data type. This should be the only C data type that needs to
be used by the end user.
* This module will almost certainly generate a compile-time warning,
such as,
subroutine lammps_open_wrapper (argc, argv, communicator, ptr) &
Variable 'communicator' at (1) is a parameter to the BIND(C)
procedure 'lammps_open_wrapper' but may not be C interoperable
This happens because lammps_open_wrapper actually takes a Fortran
INTEGER argument, whose type is defined by the MPI library itself. The
Fortran integer is converted to a C integer by the MPI library (if such
conversion is actually necessary).
* Unlike library.cpp, this module returns COPIES of the data LAMMPS actually
uses. This is done for safety reasons, as you should, in general, not be
overwriting LAMMPS data directly from Fortran. If you require this
functionality, it is possible to write another function that, for example,
returns a Fortran pointer that resolves to the C/C++ data instead of
copying the contents of that pointer to the original array as is done now.
LAMMPS.F90 defines a Fortran 2003 module, LAMMPS, which wraps all functions in
src/library.h so they can be used directly from Fortran-encoded programs.
All functions in src/library.h that use and/or return C-style pointers have
Fortran wrapper functions that use Fortran-style arrays, pointers, and
strings; all C-style memory management is handled internally with no user
intervention. See --USE-- for notes on how this interface differs from the
C interface (and the Python interface).
This interface was created by Karl Hammond who you can contact with
questions:
Karl D. Hammond
University of Tennessee, Knoxville
karlh at ugcs.caltech.edu
karlh at utk.edu
-------------------------------------
--COMPILATION--
First, be advised that mixed-language programming is not trivial. It requires
you to link in the required libraries of all languages you use (in this case,
those for Fortran, C, and C++), as well as any other libraries required.
You are also advised to read the --USE-- section below before trying to
compile.
The following steps will work to compile this module (replace ${LAMMPS_SRC}
with the path to your LAMMPS source directory).
Steps 3-5 are accomplished, possibly after some modifications to
the makefile, by make using the attached makefile. Said makefile also builds
the dynamically-linkable library (liblammps_fortran.so).
** STATIC LIBRARY INSTRUCTIONS **
(1) Compile LAMMPS as a static library.
Call the resulting file ${LAMMPS_LIB}, which will have an actual name
like liblmp_openmpi.a. If compiling using the MPI stubs in
${LAMMPS_SRC}/STUBS, you will need to know where libmpi_stubs.a
is as well (I'll call it ${MPI_STUBS} hereafter)
(2) Copy said library to your Fortran program's source directory or replace
${LAMMPS_LIB} with its full path in the instructions below.
(3) Compile (but don't link!) LAMMPS.F90. Example:
mpif90 -c LAMMPS.f90
OR
gfortran -c LAMMPS.F90
NOTE: you may get a warning such as,
subroutine lammps_open_wrapper (argc, argv, communicator, ptr) &
Variable 'communicator' at (1) is a parameter to the BIND(C)
procedure 'lammps_open_wrapper' but may not be C interoperable
This is normal (see --IMPLEMENTATION NOTES--).
(4) Compile (but don't link) LAMMPS-wrapper.cpp. You will need its header
file as well. You will have to provide the locations of LAMMPS's
header files. For example,
mpicxx -c -I${LAMMPS_SRC} LAMMPS-wrapper.cpp
OR
g++ -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp
OR
icpc -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp
(5) OPTIONAL: Make a library from the object files so you can carry around
two files instead of three. Example:
ar rs liblammps_fortran.a LAMMPS.o LAMMPS-wrapper.o
This will create the file liblammps_fortran.a that you can use in place
of "LAMMPS.o LAMMPS-wrapper.o" later. Note that you will still
need to have the .mod file from part (3).
It is also possible to add LAMMPS.o and LAMMPS-wrapper.o into the
LAMMPS library (e.g., liblmp_openmpi.a) instead of creating a separate
library, like so:
ar rs ${LAMMPS_LIB} LAMMPS.o LAMMPS-wrapper.o
In this case, you can now use the Fortran wrapper functions as if they
were part of the usual LAMMPS library interface (if you have the module
file visible to the compiler, that is).
(6) Compile (but don't link) your Fortran program. Example:
mpif90 -c myfreeformatfile.f90
mpif90 -c myfixedformatfile.f
OR
gfortran -c myfreeformatfile.f90
gfortran -c myfixedformatfile.f
The object files generated by these steps are collectively referred to
as ${my_object_files} in the next step(s).
IMPORTANT: If the Fortran module from part (3) is not in the current
directory or in one searched by the compiler for module files, you will
need to include that location via the -I flag to the compiler, like so:
mpif90 -I${LAMMPS_SRC}/examples/COUPLE/fortran2 -c myfreeformatfile.f90
(7) Link everything together, including any libraries needed by LAMMPS (such
as the C++ standard library, the C math library, the JPEG library, fftw,
etc.) For example,
mpif90 LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
${LAMMPS_LIB} -lmpi_cxx -lstdc++ -lm
OR
gfortran LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
${LAMMPS_LIB} ${MPI_STUBS} -lstdc++ -lm
OR
ifort LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
${LAMMPS_LIB} ${MPI_STUBS} -cxxlib -lm
Any other required libraries (e.g. -ljpeg, -lfftw) should be added to
the end of this line.
You should now have a working executable.
** DYNAMIC LIBRARY INSTRUCTIONS **
(1) Compile LAMMPS as a dynamic library
(make makeshlib && make -f Makefile.shlib [targetname]).
(2) Compile, but don't link, LAMMPS.F90 using the -fPIC flag, such as
mpif90 -fPIC -c LAMMPS.f90
(3) Compile, but don't link, LAMMPS-wrapper.cpp in the same manner, e.g.
mpicxx -fPIC -c LAMMPS-wrapper.cpp
(4) Make the dynamic library, like so:
mpif90 -fPIC -shared -o liblammps_fortran.so LAMMPS.o LAMMPS-wrapper.o
(5) Compile your program, such as,
mpif90 -I${LAMMPS_SRC}/examples/COUPLE/fortran2 -c myfreeformatfile.f90
where ${LAMMPS_SRC}/examples/COUPLE/fortran2 contains the .mod file from
step (3)
(6) Link everything together, such as
mpif90 ${my_object_files} -L${LAMMPS_SRC} \
-L${LAMMPS_SRC}/examples/COUPLE/fortran2 -llammps_fortran \
-llammps_openmpi -lmpi_cxx -lstdc++ -lm
If you wish to avoid the -L flags, add the directories containing your
shared libraries to the LIBRARY_PATH environment variable. At run time, you
will have to add these directories to LD_LIBRARY_PATH as well; otherwise,
your executable will not find the libraries it needs.
-------------------------------------
--USAGE--
To use this API, your program unit (PROGRAM/SUBROUTINE/FUNCTION/MODULE/etc.)
should look something like this:
program call_lammps
use LAMMPS
! Other modules, etc.
implicit none
type (lammps_instance) :: lmp ! This is a pointer to your LAMMPS instance
real (C_double) :: fix
real (C_double), dimension(:), pointer :: fix2
! Rest of declarations
call lammps_open_no_mpi ('lmp -in /dev/null -screen out.lammps',lmp)
! Set up rest of program here
call lammps_file (lmp, 'in.example')
call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1)
call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1)
call lammps_close (lmp)
end program call_lammps
Important notes:
* Though I dislike the use of pointers, they are necessary when communicating
with C and C++, which do not support Fortran's ALLOCATABLE attribute.
* There is no need to deallocate C-allocated memory; this is done for you in
the cases when it is done (which are all cases when pointers are not
accepted, such as global fix data)
* All arguments which are char* variables in library.cpp are character (len=*)
variables here. For example,
call lammps_command (lmp, 'units metal')
will work as expected.
* The public functions (the only ones you can use) have interfaces as
described in the comments at the top of LAMMPS.F90. They are not always
the same as those in library.h, since C strings are replaced by Fortran
strings and the like.
* The module attempts to check whether you have done something stupid (such
as assign a 2D array to a scalar), but it's not perfect. For example, the
command
call lammps_extract_global (nlocal, ptr, 'nlocal')
will give nlocal correctly if nlocal is a pointer to type INTEGER, but it
will give the wrong answer if nlocal is a pointer to type REAL. This is a
feature of the (void*) type cast in library.cpp. There is no way I can
check this for you! It WILL catch you if you pass it an allocatable or
fixed-size array when it expects a pointer.
* Arrays constructed from temporary data from LAMMPS are ALLOCATABLE, and
represent COPIES of data, not the originals. Functions like
lammps_extract_atom, which return actual LAMMPS data, are pointers.
* IMPORTANT: Due to the differences between C and Fortran arrays (C uses
row-major vectors, Fortran uses column-major vectors), all arrays returned
from LAMMPS have their indices swapped.
* An example of a complete program, simple.f90, is included with this
package.
-------------------------------------
--TROUBLESHOOTING--
Compile-time errors (when compiling LAMMPS.F90, that is) probably indicate
that your compiler is not new enough to support Fortran 2003 features. For
example, GCC 4.1.2 will not compile this module, but GCC 4.4.0 will.
If your compiler balks at 'use, intrinsic :: ISO_C_binding,' try removing the
intrinsic part so it looks like an ordinary module. However, it is likely
that such a compiler will also have problems with everything else in the
file as well.
If you get a segfault as soon as the lammps_open call is made, check that you
compiled your program AND LAMMPS-wrapper.cpp using the same MPI headers. Using
the stubs for one and the actual MPI library for the other will cause Bad
Things to happen.
If you find run-time errors, please pass them along via the LAMMPS Users
mailing list (please CC me as well; address above). Please provide a minimal
working example along with the names and versions of the compilers you are
using. Please make sure the error is repeatable and is in MY code, not yours
(generating a minimal working example will usually ensure this anyway).
-------------------------------------
--IMPLEMENTATION NOTES--
The Fortran procedures have the same names as the C procedures, and
their purpose is the same, but they may take different arguments. Here are
some of the important differences:
* lammps_open and lammps_open_no_mpi take a string instead of argc and
argv. This is necessary because C and C++ have a very different way
of treating strings than Fortran. If you want the command line to be
passed to lammps_open (as it often would be from C/C++), use the
GET_COMMAND intrinsic to obtain it.
* All C++ functions that accept char* pointers now accept Fortran-style
strings within this interface instead.
* All of the lammps_extract_[something] functions, which return void*
C-style pointers, have been replaced by generic subroutines that return
Fortran variables (which may be arrays). The first argument houses the
variable/pointer to be returned (pretend it's on the left-hand side); all
other arguments are identical except as stipulated above.
Note that it is not possible to declare generic functions that are selected
based solely on the type/kind/rank (TKR) signature of the return value,
only based on the TKR of the arguments.
* The SHAPE of the first argument to lammps_extract_[something] is checked
against the "shape" of the C array (e.g., double vs. double* vs. double**).
Calling a subroutine with arguments of inappropriate rank will result in an
error at run time.
* The indices i and j in lammps_extract_fix are used the same way they
are in f_ID[i][j] references in LAMMPS (i.e., starting from 1). This is
different than the way library.cpp uses these numbers, but is more
consistent with the way arrays are accessed in LAMMPS and in Fortran.
* The char* pointer normally returned by lammps_command is thrown away
in this version; note also that lammps_command is now a subroutine
instead of a function.
* The pointer to LAMMPS itself is of type(lammps_instance), which is itself
a synonym for type(C_ptr), part of ISO_C_BINDING. Type (C_ptr) is
C's void* data type.
* This module will almost certainly generate a compile-time warning,
such as,
subroutine lammps_open_wrapper (argc, argv, communicator, ptr) &
Variable 'communicator' at (1) is a parameter to the BIND(C)
procedure 'lammps_open_wrapper' but may not be C interoperable
This happens because lammps_open_wrapper actually takes a Fortran
INTEGER argument, whose type is defined by the MPI library itself. The
Fortran integer is converted to a C integer by the MPI library (if such
conversion is actually necessary).
* lammps_extract_global returns COPIES of the (scalar) data, as does the
C version.
* lammps_extract_atom, lammps_extract_compute, and lammps_extract_fix
have a first argument that will be associated with ACTUAL LAMMPS DATA.
This means the first argument must be:
* The right rank (via the DIMENSION modifier)
* A C-interoperable POINTER type (i.e., INTEGER (C_int) or
REAL (C_double)).
* lammps_extract_variable returns COPIES of the data, as the C library
interface does. There is no need to deallocate using lammps_free.
* The 'data' argument to lammps_gather_atoms and lammps_scatter atoms must
be ALLOCATABLE. It should be of type INTEGER or DOUBLE PRECISION. It
does NOT need to be C inter-operable (and indeed should not be).
* The 'count' argument of lammps_scatter_atoms is unnecessary; the shape of
the array determines the number of elements LAMMPS will read.

9
examples/COUPLE/fortran2/in.simple
View File

@ -1,10 +1,11 @@
units metal
lattice bcc 3.1656
units lj
atom_modify map array
lattice bcc 1.0
region simbox block 0 10 0 10 0 10
create_box 2 simbox
create_atoms 1 region simbox
pair_style eam/fs
pair_coeff * * path/to/my_potential.eam.fs A1 A2
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0
mass 1 58.2 # These are made-up numbers
mass 2 28.3
velocity all create 1200.0 7474848 dist gaussian

95
examples/COUPLE/fortran2/simple.f90
View File

@ -1,44 +1,111 @@
program simple
use MPI
use LAMMPS
! The following line is unnecessary, as I have included these three entities
! with the LAMMPS module, but I leave them in anyway to remind people where
! they came from
use, intrinsic :: ISO_C_binding, only : C_double, C_ptr, C_int
implicit none
type (lammps_instance) :: lmp
double precision :: compute, fix, fix2
double precision, dimension(:), allocatable :: compute_v, mass, r
double precision, dimension(:,:), allocatable :: x
real, dimension(:,:), allocatable :: x_r
! Notes:
! * If LAMMPS returns a scalar that is allocated by the library interface
! (see library.cpp), then that memory is deallocated automatically and
! the argument to lammps_extract_fix must be a SCALAR.
! * If LAMMPS returns a pointer to an array, consisting of internal LAMMPS
! data, then the argument must be an interoperable Fortran pointer.
! Interoperable means it is of type INTEGER (C_INT) or of type
! REAL (C_DOUBLE) in this context.
! * Pointers should NEVER be deallocated, as that would deallocate internal
! LAMMPS data!
! * Note that just because you can read the values of, say, a compute at
! any time does not mean those values represent the "correct" values.
! LAMMPS will abort you if you try to grab a pointer to a non-current
! entity, but once it's bound, it's your responsibility to check that
! it's current before evaluating.
! * IMPORTANT: Two-dimensional arrays (such as 'x' from extract_atom)
! will be transposed from what they might look like in C++. This is
! because of different bookkeeping conventions between Fortran and C
! that date back to about 1970 or so (when C was written).
! * Arrays start from 1, EXCEPT for mass from extract_atom, which
! starts from 0. This is because the C array actually has a blank
! first element (and thus mass[1] corresponds to the mass of type 1)
type (C_ptr) :: lmp
real (C_double), pointer :: compute => NULL()
real (C_double) :: fix, fix2
real (C_double), dimension(:), pointer :: compute_v => NULL()
real (C_double), dimension(:,:), pointer :: x => NULL()
real (C_double), dimension(:), pointer :: mass => NULL()
integer, dimension(:), allocatable :: types
double precision, dimension(:), allocatable :: r
integer :: error, narg, me, nprocs
character (len=1024) :: command_line
call MPI_Init (error)
call MPI_Comm_rank (MPI_COMM_WORLD, me, error)
call MPI_Comm_size (MPI_COMM_WORLD, nprocs, error)
! You are free to pass any string you like to lammps_open or
! lammps_open_no_mpi; here is how you pass it the command line
!call get_command (command_line)
!call lammps_open (command_line, MPI_COMM_WORLD, lmp)
! And here's how to to it with a string constant of your choice
call lammps_open_no_mpi ('lmp -log log.simple', lmp)
call lammps_open_no_mpi ('',lmp)
call lammps_file (lmp, 'in.simple')
call lammps_command (lmp, 'run 500')
! This extracts f_2 as a scalar (the last two arguments can be arbitrary)
call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1)
print *, 'Fix is ', fix
! This extracts f_4[1][1] as a scalar
call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1)
print *, 'Fix 2 is ', fix2
! This extracts the scalar compute of compute thermo_temp
call lammps_extract_compute (compute, lmp, 'thermo_temp', 0, 0)
print *, 'Compute is ', compute
! This extracts the vector compute of compute thermo_temp
call lammps_extract_compute (compute_v, lmp, 'thermo_temp', 0, 1)
print *, 'Vector is ', compute_v
! This extracts the masses
call lammps_extract_atom (mass, lmp, 'mass')
print *, 'Mass is ', mass
print *, 'Mass is ', mass(1:)
! Extracts a pointer to the arrays of positions for all atoms
call lammps_extract_atom (x, lmp, 'x')
if ( .not. allocated (x) ) print *, 'x is not allocated'
print *, 'x is ', x(1,:)
if ( .not. associated (x) ) print *, 'x is not associated'
print *, 'x is ', x(:,1) ! Prints x, y, z for atom 1
call lammps_extract_atom (x_r, lmp, 'x')
if ( .not. allocated (x_r) ) print *, 'x is not allocated'
print *, 'x_r is ', x_r(1,:)
! Extracts pointer to atom types
call lammps_gather_atoms (lmp, 'type', 1, types)
print *, 'types is ', types(1:3)
call lammps_get_coords (lmp, r)
print *, 'r is ', r(1:3)
! Allocates an array and assigns all positions to it
call lammps_gather_atoms (lmp, 'x', 3, r)
print *, 'size(r) = ', size(r)
print *, 'r is ', r(1:6)
! Puts those position data back
call lammps_scatter_atoms (lmp, 'x', r)
call lammps_command (lmp, 'run 1')
print *, 'x is ', x(:,1) ! Note that the position updates!
print *, 'Compute is ', compute ! This did only because "temp" is part of
! the thermo output; the vector part did
! not, and won't until we give LAMMPS a
! thermo output or other command that
! requires its value
call lammps_close (lmp)
call MPI_Finalize (error)
end program simple

2
python/README
View File

@ -9,7 +9,7 @@ doc/Section_python.html and in doc/Section_start.html#start_5.
Basically you need to follow these steps in the src directory:
% make makeshlib # creates Makefile.shlib
% make -f Makefile.shlib g++ # or whatever machine target you wish
% make -f Makefile.shlib g++ # build for whatever machine target you wish
% make install-python # may need to do this via sudo
You can replace the last step with running the python/install.py

19
python/examples/demo.py
View File

@ -22,11 +22,6 @@ me = 0
#nprocs = pypar.size()
from lammps import lammps
from lammps import LMPINT as INT
from lammps import LMPDOUBLE as DOUBLE
from lammps import LMPIPTR as IPTR
from lammps import LMPDPTR as DPTR
from lammps import LMPDPTRPTR as DPTRPTR
lmp = lammps()
@ -36,9 +31,9 @@ lmp.file("in.demo")
if me == 0: print "\nPython output:"
natoms = lmp.extract_global("natoms",DOUBLE)
mass = lmp.extract_atom("mass",DPTR)
x = lmp.extract_atom("x",DPTRPTR)
natoms = lmp.extract_global("natoms",0)
mass = lmp.extract_atom("mass",2)
x = lmp.extract_atom("x",3)
print "Natoms, mass, x[0][0] coord =",natoms,mass[1],x[0][0]
temp = lmp.extract_compute("thermo_temp",0,0)
@ -53,13 +48,13 @@ print "Velocity component from atom-style variable =",vy[1]
natoms = lmp.get_natoms()
print "Natoms from get_natoms =",natoms
xc = lmp.get_coords()
print "Global coords from get_coords =",xc[0],xc[1],xc[31]
xc = lmp.gather_atoms("x",1,3)
print "Global coords from gather_atoms =",xc[0],xc[1],xc[31]
xc[0] = xc[0] + 1.0
lmp.put_coords(xc)
lmp.scatter_atoms("x",1,3,xc)
print "Changed x[0][0] via put_coords =",x[0][0]
print "Changed x[0][0] via scatter_atoms =",x[0][0]
# uncomment if running in parallel via Pypar
#print "Proc %d out of %d procs has" % (me,nprocs), lmp

12
python/lammps.py
View File

@ -24,8 +24,8 @@ class lammps:
# if name = "g++", load liblammps_g++.so
try:
if not name: self.lib = CDLL("liblammps.so")
else: self.lib = CDLL("liblammps_%s.so" % name)
if not name: self.lib = CDLL("liblammps.so",RTLD_GLOBAL)
else: self.lib = CDLL("liblammps_%s.so" % name,RTLD_GLOBAL)
except:
type,value,tb = sys.exc_info()
traceback.print_exception(type,value,tb)
@ -89,11 +89,11 @@ class lammps:
self.lib.lammps_extract_compute.restype = POINTER(c_double)
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
return ptr[0]
elif type == 1:
if type == 1:
self.lib.lammps_extract_compute.restype = POINTER(c_double)
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
return ptr
elif type == 2:
if type == 2:
self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double))
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
return ptr
@ -110,11 +110,11 @@ class lammps:
result = ptr[0]
self.lib.lammps_free(ptr)
return result
elif type == 1:
if type == 1:
self.lib.lammps_extract_fix.restype = POINTER(c_double)
ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j)
return ptr
elif type == 2:
if type == 2:
self.lib.lammps_extract_fix.restype = POINTER(POINTER(c_double))
ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j)
return ptr

4
src/Makefile.shlib
View File

File diff suppressed because one or more lines are too long

12
src/REPLICA/verlet_split.cpp
View File

@ -1,4 +1,4 @@
/* -------------------------------------------------------------------------
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -274,6 +274,13 @@ void VerletSplit::run(int n)
rk_setup();
// check if OpenMP support fix defined
Fix *fix_omp;
int ifix = modify->find_fix("package_omp");
if (ifix < 0) fix_omp = NULL;
else fix_omp = modify->fix[ifix];
// flags for timestepping iterations
int n_post_integrate = modify->n_post_integrate;
@ -372,7 +379,8 @@ void VerletSplit::run(int n)
} else {
// run FixOMP as sole pre_force fix, if defined
if (fix_omp != NULL) fix_omp->pre_force(vflag);
if (fix_omp) fix_omp->pre_force(vflag);
if (force->kspace) {
timer->stamp();

64
src/RIGID/fix_rigid.cpp
View File

@ -695,39 +695,20 @@ void FixRigid::setup(int vflag)
tagint *image = atom->image;
double **x = atom->x;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double xy = domain->xy;
double xz = domain->xz;
double yz = domain->yz;
double dx,dy,dz;
double unwrap[3];
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
int xbox,ybox,zbox;
double xunwrap,yunwrap,zunwrap,dx,dy,dz;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
} else {
xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
yunwrap = x[i][1] + ybox*yprd + zbox*yz;
zunwrap = x[i][2] + zbox*zprd;
}
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
@ -929,17 +910,9 @@ void FixRigid::final_integrate()
double **f = atom->f;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
if (triclinic) {
xy = domain->xy;
xz = domain->xz;
yz = domain->yz;
}
double dx,dy,dz;
double unwrap[3];
int xbox,ybox,zbox;
double xunwrap,yunwrap,zunwrap,dx,dy,dz;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
@ -951,23 +924,10 @@ void FixRigid::final_integrate()
sum[ibody][1] += f[i][1];
sum[ibody][2] += f[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
} else {
xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
yunwrap = x[i][1] + ybox*yprd + zbox*yz;
zunwrap = x[i][2] + zbox*zprd;
}
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2];
sum[ibody][3] += dy*f[i][2] - dz*f[i][1];
sum[ibody][4] += dz*f[i][0] - dx*f[i][2];

42
src/USER-MISC/compute_temp_rotate.cpp
View File

@ -89,11 +89,8 @@ double ComputeTempRotate::compute_scalar()
{
double vthermal[3];
double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3];
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
invoked_scalar = update->ntimestep;
@ -123,17 +120,13 @@ double ComputeTempRotate::compute_scalar()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2];
vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0];
vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1];
vthermal[0] = v[i][0] - vbiasall[i][0];
vthermal[1] = v[i][1] - vbiasall[i][1];
vthermal[2] = v[i][2] - vbiasall[i][2];
@ -155,14 +148,10 @@ double ComputeTempRotate::compute_scalar()
void ComputeTempRotate::compute_vector()
{
int i;
double vthermal[3];
double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3];
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
invoked_vector = update->ntimestep;
@ -189,25 +178,20 @@ void ComputeTempRotate::compute_vector()
}
double massone,t[6];
for (i = 0; i < 6; i++) t[i] = 0.0;
for (int i = 0; i < 6; i++) t[i] = 0.0;
for (i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2];
vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0];
vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1];
vthermal[0] = v[i][0] - vbiasall[i][0];
vthermal[1] = v[i][1] - vbiasall[i][1];
vthermal[2] = v[i][2] - vbiasall[i][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0];
@ -219,7 +203,7 @@ void ComputeTempRotate::compute_vector()
}
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ----------------------------------------------------------------------

37
src/USER-MISC/fix_addtorque.cpp
View File

@ -109,19 +109,22 @@ void FixAddTorque::init()
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist");
if (xvar < 0)
error->all(FLERR,"Variable name for fix addtorque does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else error->all(FLERR,"Variable for fix addtorque is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist");
if (yvar < 0)
error->all(FLERR,"Variable name for fix addtorque does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else error->all(FLERR,"Variable for fix addtorque is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all(FLERR,"Variable name for fix addtorque does not exist");
if (zvar < 0)
error->all(FLERR,"Variable name for fix addtorque does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else error->all(FLERR,"Variable for fix addtorque is invalid style");
}
@ -168,16 +171,14 @@ void FixAddTorque::post_force(int vflag)
int nlocal = atom->nlocal;
double mvv2e = force->mvv2e;
int xbox,ybox,zbox;
double dx,dy,dz,vx,vy,vz,fx,fy,fz,massone,omegadotr;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double tcm[3],xcm[3],angmom[3],omega[3],itorque[3],domegadt[3],tlocal[3];
double inertia[3][3];
double unwrap[3];
// foriginal[0] = "potential energy" for added force
// foriginal[123] = torque on atoms before extra force added
foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
force_flag = 0;
@ -200,12 +201,10 @@ void FixAddTorque::post_force(int vflag)
tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
omegadotr = omega[0]*dx+omega[1]*dy+omega[2]*dz;
@ -222,12 +221,10 @@ void FixAddTorque::post_force(int vflag)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
vx = mvv2e*(dz*omega[1]-dy*omega[2]);
vy = mvv2e*(dx*omega[2]-dz*omega[0]);
vz = mvv2e*(dy*omega[0]-dx*omega[1]);

16
src/compute_com_molecule.cpp
View File

@ -96,8 +96,8 @@ void ComputeCOMMolecule::init()
void ComputeCOMMolecule::compute_array()
{
int i,imol;
double xbox,ybox,zbox;
double massone;
double unwrap[3];
invoked_array = update->ntimestep;
@ -113,23 +113,17 @@ void ComputeCOMMolecule::compute_array()
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
com[imol][0] += (x[i][0] + xbox*xprd) * massone;
com[imol][1] += (x[i][1] + ybox*yprd) * massone;
com[imol][2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,

17
src/compute_gyration.cpp
View File

@ -83,25 +83,22 @@ void ComputeGyration::compute_vector()
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double rg[6];
rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
rg[0] += dx*dx * massone;
rg[1] += dy*dy * massone;
rg[2] += dz*dz * massone;

57
src/compute_gyration_molecule.cpp
View File

@ -123,8 +123,8 @@ void ComputeGyrationMolecule::init()
void ComputeGyrationMolecule::compute_vector()
{
int i,imol;
double xbox,ybox,zbox,dx,dy,dz;
double massone;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
@ -141,21 +141,16 @@ void ComputeGyrationMolecule::compute_vector()
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
dx = (x[i][0] + xbox*xprd) - comall[imol][0];
dy = (x[i][1] + ybox*yprd) - comall[imol][1];
dz = (x[i][2] + zbox*zprd) - comall[imol][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[imol][0];
dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rg[imol] += (dx*dx + dy*dy + dz*dz) * massone;
@ -171,8 +166,8 @@ void ComputeGyrationMolecule::compute_vector()
void ComputeGyrationMolecule::compute_array()
{
int i,j,imol;
double xbox,ybox,zbox,dx,dy,dz;
double massone;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
@ -191,21 +186,15 @@ void ComputeGyrationMolecule::compute_array()
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
dx = (x[i][0] + xbox*xprd) - comall[imol][0];
dy = (x[i][1] + ybox*yprd) - comall[imol][1];
dz = (x[i][2] + zbox*zprd) - comall[imol][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - comall[imol][0];
dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rgt[imol][0] += dx*dx * massone;
@ -233,8 +222,8 @@ void ComputeGyrationMolecule::compute_array()
void ComputeGyrationMolecule::molcom()
{
int i,imol;
double xbox,ybox,zbox,dx,dy,dz;
double massone;
double dx,dy,dz,massone;
double unwrap[3];
for (i = 0; i < nmolecules; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0;
@ -248,23 +237,17 @@ void ComputeGyrationMolecule::molcom()
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
com[imol][0] += (x[i][0] + xbox*xprd) * massone;
com[imol][1] += (x[i][1] + ybox*yprd) * massone;
com[imol][2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,

22
src/compute_msd_molecule.cpp
View File

@ -111,8 +111,8 @@ void ComputeMSDMolecule::init()
void ComputeMSDMolecule::compute_array()
{
int i,imol;
double xbox,ybox,zbox,dx,dy,dz;
double massone;
double dx,dy,dz,massone;
double unwrap[3];
invoked_array = update->ntimestep;
@ -130,23 +130,17 @@ void ComputeMSDMolecule::compute_array()
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i];
if (molmap) imol = molmap[imol-idlo];
else imol--;
com[imol][0] += (x[i][0] + xbox*xprd) * massone;
com[imol][1] += (x[i][1] + ybox*yprd) * massone;
com[imol][2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
}
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,

15
src/fix_addforce.cpp
View File

@ -217,12 +217,6 @@ void FixAddForce::min_setup(int vflag)
void FixAddForce::post_force(int vflag)
{
int xbox,ybox,zbox;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
@ -247,18 +241,15 @@ void FixAddForce::post_force(int vflag)
// potential energy = - x dot f in unwrapped coords
if (varflag == CONSTANT) {
double unwrap[3];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (iregion >= 0 &&
!domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
continue;
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
foriginal[0] -= xvalue * (x[i][0]+xbox*xprd) +
yvalue * (x[i][1]+ybox*yprd) + zvalue * (x[i][2]+zbox*zprd);
domain->unmap(x[i],image[i],unwrap);
foriginal[0] -= xvalue*unwrap[0] + yvalue*unwrap[1] + zvalue*unwrap[2];
foriginal[1] += f[i][0];
foriginal[2] += f[i][1];
foriginal[3] += f[i][2];

18
src/fix_momentum.cpp
View File

@ -57,7 +57,8 @@ FixMomentum::FixMomentum(LAMMPS *lmp, int narg, char **arg) :
if (linear)
if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 ||
zflag < 0 || zflag > 1) error->all(FLERR,"Illegal fix momentum command");
zflag < 0 || zflag > 1)
error->all(FLERR,"Illegal fix momentum command");
// cannot have 0 atoms in group
@ -121,20 +122,15 @@ void FixMomentum::end_of_step()
tagint *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
v[i][0] -= omega[1]*dz - omega[2]*dy;
v[i][1] -= omega[2]*dx - omega[0]*dz;
v[i][2] -= omega[0]*dy - omega[1]*dx;

15
src/fix_spring_rg.cpp
View File

@ -110,20 +110,15 @@ void FixSpringRG::post_force(int vflag)
int nlocal = atom->nlocal;
double massfrac;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
int xbox,ybox,zbox;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
term1 = 2.0 * k * (1.0 - rg0/rg);
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
massfrac = mass[type[i]]/masstotal;
f[i][0] -= term1*dx*massfrac;
f[i][1] -= term1*dy*massfrac;

34
src/fix_spring_self.cpp
View File

@ -45,9 +45,10 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command");
xflag = yflag = zflag = 1;
if (narg == 5) {
if (strcmp(arg[4],"xyz") == 0) {
; /* default */
xflag = yflag = zflag = 1;
} else if (strcmp(arg[4],"xy") == 0) {
zflag = 0;
} else if (strcmp(arg[4],"xz") == 0) {
@ -78,20 +79,9 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
tagint *image = atom->image;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xoriginal[i][0] = x[i][0] + xbox*xprd;
xoriginal[i][1] = x[i][1] + ybox*yprd;
xoriginal[i][2] = x[i][2] + zbox*zprd;
} else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
}
espring = 0.0;
@ -161,21 +151,17 @@ void FixSpringSelf::post_force(int vflag)
tagint *image = atom->image;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
double dx,dy,dz;
double unwrap[3];
espring = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = x[i][0] + xbox*xprd - xoriginal[i][0];
dy = x[i][1] + ybox*yprd - xoriginal[i][1];
dz = x[i][2] + zbox*zprd - xoriginal[i][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xoriginal[i][0];
dy = unwrap[1] - xoriginal[i][1];
dz = unwrap[2] - xoriginal[i][2];
if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0;

51
src/fix_tmd.cpp
View File

@ -98,26 +98,17 @@ FixTMD::FixTMD(LAMMPS *lmp, int narg, char **arg) :
double *mass = atom->mass;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double dx,dy,dz;
int xbox,ybox,zbox;
rho_start = 0.0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = x[i][0] + xbox*xprd - xf[i][0];
dy = x[i][1] + ybox*yprd - xf[i][1];
dz = x[i][2] + zbox*zprd - xf[i][2];
domain->unmap(x[i],image[i],xold[i]);
dx = xold[i][0] - xf[i][0];
dy = xold[i][1] - xf[i][1];
dz = xold[i][2] - xf[i][2];
rho_start += mass[type[i]]*(dx*dx + dy*dy + dz*dz);
xold[i][0] = x[i][0] + xbox*xprd;
xold[i][1] = x[i][1] + ybox*yprd;
xold[i][2] = x[i][2] + zbox*zprd;
} else xold[i][0] = xold[i][1] = xold[i][2] = 0.0;
}
@ -190,10 +181,8 @@ void FixTMD::initial_integrate(int vflag)
double dxold,dyold,dzold,xback,yback,zback;
double gamma_forward,gamma_back,gamma_max,lambda;
double kt,fr,kttotal,frtotal,dtfm;
double unwrap[3];
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@ -202,7 +191,6 @@ void FixTMD::initial_integrate(int vflag)
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
@ -216,12 +204,10 @@ void FixTMD::initial_integrate(int vflag)
dxold = xold[i][0] - xf[i][0];
dyold = xold[i][1] - xf[i][1];
dzold = xold[i][2] - xf[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = x[i][0] + xbox*xprd - xf[i][0];
dy = x[i][1] + ybox*yprd - xf[i][1];
dz = x[i][2] + zbox*zprd - xf[i][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xf[i][0];
dy = unwrap[1] - xf[i][1];
dz = unwrap[2] - xf[i][2];
a += mass[type[i]]*(dxold*dxold + dyold*dyold + dzold*dzold);
b += mass[type[i]]*(dx *dxold + dy *dyold + dz *dzold);
e += mass[type[i]]*(dx *dx + dy *dy + dz *dz);
@ -258,12 +244,10 @@ void FixTMD::initial_integrate(int vflag)
dxold = xold[i][0] - xf[i][0];
dyold = xold[i][1] - xf[i][1];
dzold = xold[i][2] - xf[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xback = x[i][0] + xbox*xprd + gamma_back*dxold;
yback = x[i][1] + ybox*yprd + gamma_back*dyold;
zback = x[i][2] + zbox*zprd + gamma_back*dzold;
domain->unmap(x[i],image[i],unwrap);
xback = unwrap[0] + gamma_back*dxold;
yback = unwrap[1] + gamma_back*dyold;
zback = unwrap[2] + gamma_back*dzold;
dxkt = xback - xold[i][0];
dykt = yback - xold[i][1];
dzkt = zback - xold[i][2];
@ -314,12 +298,7 @@ void FixTMD::initial_integrate(int vflag)
x[i][2] += gamma_forward*dzold;
v[i][2] += gamma_forward*dzold/dtv;
f[i][2] += gamma_forward*dzold/dtv/dtfm;
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xold[i][0] = x[i][0] + xbox*xprd;
xold[i][1] = x[i][1] + ybox*yprd;
xold[i][2] = x[i][2] + zbox*zprd;
domain->unmap(x[i],image[i],xold[i]);
}
}
}

176
src/group.cpp
View File

@ -766,34 +766,27 @@ void Group::xcm(int igroup, double masstotal, double *cm)
double cmone[3];
cmone[0] = cmone[1] = cmone[2] = 0.0;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
double massone;
double unwrap[3];
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = rmass[i];
cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
cmone[0] += unwrap[0] * massone;
cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = mass[type[i]];
cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
cmone[0] += unwrap[0] * massone;
cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
}
}
@ -827,34 +820,27 @@ void Group::xcm(int igroup, double masstotal, double *cm, int iregion)
double cmone[3];
cmone[0] = cmone[1] = cmone[2] = 0.0;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
double massone;
double unwrap[3];
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = rmass[i];
cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
cmone[0] += unwrap[0] * massone;
cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = mass[type[i]];
cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone;
domain->unmap(x[i],image[i],unwrap);
cmone[0] += unwrap[0] * massone;
cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
}
}
@ -1102,21 +1088,16 @@ double Group::gyration(int igroup, double masstotal, double *cm)
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double rg = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rg += (dx*dx + dy*dy + dz*dz) * massone;
@ -1147,21 +1128,16 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion)
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double rg = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rg += (dx*dx + dy*dy + dz*dz) * massone;
@ -1192,22 +1168,18 @@ void Group::angmom(int igroup, double *cm, double *lmom)
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double p[3];
p[0] = p[1] = p[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
@ -1238,22 +1210,18 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double p[3];
p[0] = p[1] = p[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
@ -1280,22 +1248,18 @@ void Group::torque(int igroup, double *cm, double *tq)
tagint *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0];
@ -1321,22 +1285,18 @@ void Group::torque(int igroup, double *cm, double *tq, int iregion)
tagint *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0];
@ -1364,11 +1324,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double ione[3][3];
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
@ -1376,12 +1334,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
ione[0][0] += massone * (dy*dy + dz*dz);
@ -1418,11 +1374,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
double ione[3][3];
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
@ -1430,12 +1384,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - cm[0];
dy = unwrap[1] - cm[1];
dz = unwrap[2] - cm[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
ione[0][0] += massone * (dy*dy + dz*dz);

15
src/velocity.cpp
View File

@ -718,20 +718,15 @@ void Velocity::zero_rotation()
tagint *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double unwrap[3];
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0]- xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
v[i][0] -= omega[1]*dz - omega[2]*dy;
v[i][1] -= omega[2]*dx - omega[0]*dz;
v[i][2] -= omega[0]*dy - omega[1]*dx;

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "23 Oct 2012"
#define LAMMPS_VERSION "26 Oct 2012"