Merge branch 'master' into progguide-axel

This commit is contained in:
Axel Kohlmeyer
2020-10-08 21:51:22 -04:00
18 changed files with 14469 additions and 1221 deletions

View File

@ -81,6 +81,7 @@ OPT.
* :doc:`coul/slater/long <pair_coul_slater>`
* :doc:`coul/shield <pair_coul_shield>`
* :doc:`coul/streitz <pair_coul>`
* :doc:`coul/tt <pair_coul_tt>`
* :doc:`coul/wolf (ko) <pair_coul>`
* :doc:`coul/wolf/cs <pair_cs>`
* :doc:`dpd (gio) <pair_dpd>`

View File

@ -24,7 +24,7 @@ of the source files: the lammps.f90 file needs to be compiled first,
since it provides the ``LIBLAMMPS`` module that is imported by the
Fortran code using the interface.
.. versionadded:: 30Sep2020
.. versionadded:: 6Oct2020
.. admonition:: Work in Progress
:class: note

131
doc/src/pair_coul_tt.rst Normal file
View File

@ -0,0 +1,131 @@
.. index:: pair_style coul/tt
pair_style coul/tt command
==========================
Syntax
""""""
.. code-block:: LAMMPS
pair_style style args
* style = *coul/tt*
* args = list of arguments for a particular style
.. parsed-literal::
*coul/tt* args = n cutoff
n = degree of polynomial
cutoff = global cutoff (distance units)
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay ... coul/tt 4 12.0
pair_coeff 1 2 coul/tt 4.5 1.0
pair_coeff 1 2 coul/tt 4.0 1.0 4 12.0
pair_coeff 1 3* coul/tt 4.5 1.0 4
Description
"""""""""""
The *coul/tt* pair style is meant to be used with force fields that
include explicit polarization through Drude dipoles.
The *coul/tt* pair style should be used as a sub-style within in the
:doc:`pair_style hybrid/overlay <pair_hybrid>` command, in conjunction with a
main pair style including Coulomb interactions and *thole* pair style,
or with *lj/cut/thole/long* pair style that is equivalent to the combination
of preceding two.
The *coul/tt* pair styles compute the charge-dipole Coulomb interaction damped
at short distances by a function
.. math::
f_{n,ij}(r) = 1 - c_{ij} \cdot e^{-b_{ij} r} \sum_{k=0}^n \frac{(b_{ij} r)^k}{k!}
This function results from an adaptation to the Coulomb interaction :ref:`(Salanne)
<Salanne1>` of the damping function originally proposed
by :ref:`Tang Toennies <TangToennies1>` for van der Waals interactions.
The polynomial takes the degree 4 for damping the Coulomb interaction.
The parameters :math:`b_{ij}` and :math:`c_{ij}` could be determined from
first-principle calculations for small, mainly mono-atomic, ions :ref:`(Salanne)
<Salanne1>`, or else treated as empirical for large molecules.
In pair styles with Drude induced dipoles, this damping function is typically
applied to the interactions between a Drude charge (either :math:`q_{D,i}` on
a Drude particle or :math:`-q_{D,i}` on the respective
Drude core)) and a charge on a non-polarizable atom, :math:`q_{j}`.
The Tang-Toennies function could also be used to damp electrostatic
interactions between the (non-polarizable part of the) charge of a core,
:math:`q_{i}-q_{D,i}`, and the Drude charge of another, :math:`-q_{D,j}`.
The :math:`b_{ij}` and :math:`c_{ij}` are equal to :math:`b_{ji}` and
:math:`c_{ji}` in the case of core-core interactions.
For pair_style *coul/tt*\ , the following coefficients must be defined for
each pair of atoms types via the :doc:`pair_coeff <pair_coeff>` command
as in the example above.
* :math:`b_{ij}`
* :math:`c_{ij}`
* degree of polynomial (positive integer)
* cutoff (distance units)
The last two coefficients are optional. If not specified the global
degree of the polynomial or the global cutoff specified in the pair_style
command are used. In order to specify a cutoff (forth argument), the degree of
the polynomial (third argument) must also be specified.
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
The *coul/tt* pair style does not support mixing. Thus, coefficients
for all I,J pairs must be specified explicitly.
Restrictions
""""""""""""
These pair styles are part of the USER-DRUDE package. They are only
enabled if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` doc page for more info.
This pair_style should currently not be used with the :doc:`charmm dihedral
style <dihedral_charmm>` if the latter has non-zero 1-4 weighting
factors. This is because the *coul/tt* pair style does not know which
pairs are 1-4 partners of which dihedrals.
Related commands
""""""""""""""""
:doc:`fix drude <fix_drude>`, :doc:`fix langevin/drude <fix_langevin_drude>`,
:doc:`fix drude/transform <fix_drude_transform>`,
:doc:`compute temp/drude <compute_temp_drude>`,
:doc:`pair_style thole <pair_thole>`
Default
"""""""
none
----------
.. _Thole1:
**(Thole)** Chem Phys, 59, 341 (1981).
.. _Salanne1:
**(Salanne)** Salanne, Rotenberg, Jahn, Vuilleumier, Simon, Christian and Madden, Theor Chem Acc, 131, 1143 (2012).
.. _TangToennies1:
**(Tang and Toennies)** J Chem Phys, 80, 3726 (1984).

View File

@ -144,6 +144,7 @@ accelerated styles exist.
* :doc:`coul/slater/long <pair_coul>` - long-range smeared out Coulomb
* :doc:`coul/shield <pair_coul_shield>` - Coulomb for boron nitride for use with :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` potential
* :doc:`coul/streitz <pair_coul>` - Coulomb via Streitz/Mintmire Slater orbitals
* :doc:`coul/tt <pair_coul_tt>` - damped charge-dipole Coulomb for Drude dipoles
* :doc:`coul/wolf <pair_coul>` - Coulomb via Wolf potential
* :doc:`coul/wolf/cs <pair_cs>` - Coulomb via Wolf potential with core/shell adjustments
* :doc:`dpd <pair_dpd>` - dissipative particle dynamics (DPD)

View File

@ -1410,6 +1410,7 @@ Izz
Jacobsen
Jadhav
jagreat
Jahn
Jalalvand
james
Janssen
@ -2715,6 +2716,7 @@ Rosenberger
Rossky
rosybrown
rotationally
Rotenberg
Rovigatti
royalblue
rozero
@ -2749,6 +2751,7 @@ safezone
Safran
Sagui
Saidi
Salanne
Salles
sandia
Sandia
@ -3112,6 +3115,7 @@ Tohoku
tokenizer
tokyo
tol
Toennies
tomic
toolchain
topologies
@ -3340,6 +3344,7 @@ vtp
vtr
vtu
Vu
Vuilleumier
vv
vx
Vx

File diff suppressed because it is too large Load Diff

View File

@ -7,8 +7,9 @@ LAMMPS_SRC = $(LAMMPS_ROOT)/src
# Uncomment the line below if using the MPI stubs library
MPI_STUBS = #-I$(LAMMPS_SRC)/STUBS
FC = mpif90 # replace with your Fortran compiler
CXX = mpicxx # replace with your C++ compiler
FC = mpifort # replace with your Fortran compiler
CXX = mpicxx # replace with your C++ compiler
CXXLIB = -lstdc++ # replace with your C++ runtime libs
# Flags for Fortran compiler, C++ compiler, and C preprocessor, respectively
FFLAGS = -O2 -fPIC
@ -18,7 +19,7 @@ CPPFLAGS = -DOMPI_SKIP_MPICXX=1 -DMPICH_SKIP_MPICXX
all : liblammps_fortran.a liblammps_fortran.so
liblammps_fortran.so : LAMMPS.o LAMMPS-wrapper.o
$(FC) $(FFLAGS) -shared -o $@ $^
$(FC) $(FFLAGS) -shared -o $@ $^ $(CXXLIB)
liblammps_fortran.a : LAMMPS.o LAMMPS-wrapper.o
$(AR) rs $@ $^

View File

@ -1,265 +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. 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.
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.

View File

@ -60,19 +60,19 @@ program 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)
call lammps_extract_fix (fix, lmp, '2', LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR, 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)
call lammps_extract_fix (fix2, lmp, '4', LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY, 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)
call lammps_extract_compute (compute, lmp, 'thermo_temp', LMP_STYLE_GLOBAL, LMP_TYPE_SCALAR)
print *, 'Compute is ', compute
! This extracts the vector compute of compute thermo_temp
call lammps_extract_compute (compute_v, lmp, 'thermo_temp', 0, 1)
call lammps_extract_compute (compute_v, lmp, 'thermo_temp', LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR)
print *, 'Vector is ', compute_v
! This extracts the masses
@ -90,6 +90,7 @@ program simple
! Allocates an array and assigns all positions to it
call lammps_gather_atoms (lmp, 'x', 3, r)
print *, 'natoms = ', int(lammps_get_natoms(lmp))
print *, 'size(r) = ', size(r)
print *, 'r is ', r(1:6)

View File

@ -14,3 +14,6 @@ use of the `extra/special/per/atom` keyword
* `swm4-ndp` -- 4-site rigid water model in NpT ensemble (no Thole
damping)
* `ethylene_glycol` -- simulation in NpT ensemble with Thole
and Tang-Toennies damping and with Nosé-Hoover thermostat

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,80 @@
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
pair_style hybrid/overlay lj/cut/coul/long 8.0 8.0 thole 2.600 8.0 coul/tt 4 8.0
pair_modify tail yes
kspace_style pppm 1.0e-5
read_data data.ethylene_glycol
pair_coeff 1 1 lj/cut/coul/long 0.057289 3.500000 # CTO CTO ~
pair_coeff 1 2 lj/cut/coul/long 0.091945 3.304542 # CTO OHG ~
pair_coeff 1 3 lj/cut/coul/long 0.038625 2.958040 # CTO H1O ~
pair_coeff 1 4 lj/cut/coul/long 0.000000 0.000000 # CTO HOG ~
pair_coeff 2 2 lj/cut/coul/long 0.147565 3.120000 # OHG OHG ~
pair_coeff 2 3 lj/cut/coul/long 0.061990 2.792848 # OHG H1O ~
pair_coeff 2 4 lj/cut/coul/long 0.000000 0.000000 # OHG HOG ~
pair_coeff 3 3 lj/cut/coul/long 0.026041 2.500000 # H1O H1O ~
pair_coeff 3 4 lj/cut/coul/long 0.000000 0.000000 # H1O HOG ~
pair_coeff 4 4 lj/cut/coul/long 0.000000 0.000000 # HOG HOG ~
pair_coeff * 5* lj/cut/coul/long 0.000000 0.000000
pair_coeff 1 1 thole 1.662
pair_coeff 1 2 thole 1.561
pair_coeff 1 5 thole 1.662
pair_coeff 1 6 thole 1.561
pair_coeff 2 2 thole 1.467
pair_coeff 2 5 thole 1.561
pair_coeff 2 6 thole 1.467
pair_coeff 5 5 thole 1.662
pair_coeff 5 6 thole 1.561
pair_coeff 6 6 thole 1.467
pair_coeff 2 4 coul/tt 4.5 1.0
pair_coeff 4 6 coul/tt 4.5 1.0
pair_coeff 1 4 coul/tt 4.5 1.0
pair_coeff 4 5 coul/tt 4.5 1.0
group ATOMS type 1 2 3 4
group CORES type 1 2
group DRUDES type 5 6
fix DRUDE all drude C C N N D D
fix SHAKE ATOMS shake 0.0001 20 0 b 3 4
neighbor 2.0 bin
timestep 1.0
variable TK equal 298.0
variable TDK equal 1.0
variable PBAR equal 1.0
comm_modify vel yes
velocity ATOMS create ${TK} 12345
compute TATOM ATOMS temp
compute TDRUDE all temp/drude
fix DTDIR all drude/transform/direct
fix TSTAT ATOMS npt temp ${TK} ${TK} 200 iso ${PBAR} ${PBAR} 1000
fix_modify TSTAT temp TATOM press thermo_press
fix TSTDR DRUDES nvt temp ${TDK} ${TDK} 50
fix DTINV all drude/transform/inverse
fix ICECUBE all momentum 1000 linear 1 1 1
thermo_style custom step etotal ke pe ebond eangle evdwl ecoul elong &
press vol density c_TATOM c_TDRUDE[1] c_TDRUDE[2]
thermo 50
run 500

View File

@ -0,0 +1,200 @@
LAMMPS (24 Aug 2020)
using 1 OpenMP thread(s) per MPI task
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
pair_style hybrid/overlay lj/cut/coul/long 8.0 8.0 thole 2.600 8.0 coul/tt 4 8.0
pair_modify tail yes
kspace_style pppm 1.0e-5
read_data data.ethylene_glycol
Reading data file ...
orthogonal box = (0.0000000 0.0000000 0.0000000) to (35.000000 35.000000 35.000000)
1 by 1 by 1 MPI processor grid
reading atoms ...
2800 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
9 = max dihedrals/atom
reading bonds ...
2600 bonds
reading angles ...
2800 angles
reading dihedrals ...
3000 dihedrals
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0.0 0.0 0.5
special bond factors coul: 0.0 0.0 0.5
5 = max # of 1-2 neighbors
6 = max # of 1-3 neighbors
10 = max # of 1-4 neighbors
13 = max # of special neighbors
special bonds CPU = 0.002 seconds
read_data CPU = 0.023 seconds
pair_coeff 1 1 lj/cut/coul/long 0.057289 3.500000 # CTO CTO ~
pair_coeff 1 2 lj/cut/coul/long 0.091945 3.304542 # CTO OHG ~
pair_coeff 1 3 lj/cut/coul/long 0.038625 2.958040 # CTO H1O ~
pair_coeff 1 4 lj/cut/coul/long 0.000000 0.000000 # CTO HOG ~
pair_coeff 2 2 lj/cut/coul/long 0.147565 3.120000 # OHG OHG ~
pair_coeff 2 3 lj/cut/coul/long 0.061990 2.792848 # OHG H1O ~
pair_coeff 2 4 lj/cut/coul/long 0.000000 0.000000 # OHG HOG ~
pair_coeff 3 3 lj/cut/coul/long 0.026041 2.500000 # H1O H1O ~
pair_coeff 3 4 lj/cut/coul/long 0.000000 0.000000 # H1O HOG ~
pair_coeff 4 4 lj/cut/coul/long 0.000000 0.000000 # HOG HOG ~
pair_coeff * 5* lj/cut/coul/long 0.000000 0.000000
pair_coeff 1 1 thole 1.662
pair_coeff 1 2 thole 1.561
pair_coeff 1 5 thole 1.662
pair_coeff 1 6 thole 1.561
pair_coeff 2 2 thole 1.467
pair_coeff 2 5 thole 1.561
pair_coeff 2 6 thole 1.467
pair_coeff 5 5 thole 1.662
pair_coeff 5 6 thole 1.561
pair_coeff 6 6 thole 1.467
pair_coeff 2 4 coul/tt 4.5 1.0
pair_coeff 4 6 coul/tt 4.5 1.0
pair_coeff 1 4 coul/tt 4.5 1.0
pair_coeff 4 5 coul/tt 4.5 1.0
group ATOMS type 1 2 3 4
2000 atoms in group ATOMS
group CORES type 1 2
800 atoms in group CORES
group DRUDES type 5 6
800 atoms in group DRUDES
fix DRUDE all drude C C N N D D
fix SHAKE ATOMS shake 0.0001 20 0 b 3 4
400 = # of size 2 clusters
400 = # of size 3 clusters
0 = # of size 4 clusters
0 = # of frozen angles
find clusters CPU = 0.001 seconds
neighbor 2.0 bin
timestep 1.0
variable TK equal 298.0
variable TDK equal 1.0
variable PBAR equal 1.0
comm_modify vel yes
velocity ATOMS create ${TK} 12345
velocity ATOMS create 298 12345
compute TATOM ATOMS temp
compute TDRUDE all temp/drude
fix DTDIR all drude/transform/direct
fix TSTAT ATOMS npt temp ${TK} ${TK} 200 iso ${PBAR} ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 ${TK} 200 iso ${PBAR} ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 298 200 iso ${PBAR} ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 298 200 iso 1 ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 298 200 iso 1 1 1000
fix_modify TSTAT temp TATOM press thermo_press
WARNING: Temperature for fix modify is not for group all (src/fix_nh.cpp:1428)
fix TSTDR DRUDES nvt temp ${TDK} ${TDK} 50
fix TSTDR DRUDES nvt temp 1 ${TDK} 50
fix TSTDR DRUDES nvt temp 1 1 50
fix DTINV all drude/transform/inverse
fix ICECUBE all momentum 1000 linear 1 1 1
thermo_style custom step etotal ke pe ebond eangle evdwl ecoul elong press vol density c_TATOM c_TDRUDE[1] c_TDRUDE[2]
thermo 50
run 500
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:328)
G vector (1/distance) = 0.41206781
grid = 54 54 54
stencil order = 5
estimated absolute RMS force accuracy = 0.0040479865
estimated relative force accuracy = 1.2190391e-05
using double precision FFTW3
3d grid and FFT values/proc = 250047 157464
Rebuild special list taking Drude particles into account
Old max number of 1-2 to 1-4 neighbors: 13
New max number of 1-2 to 1-4 neighbors: 13 (+0)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10
ghost atom cutoff = 10
binsize = 5, bins = 7 7 7
3 neighbor lists, perpetual/occasional/extra = 3 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair thole, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(3) pair coul/tt, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 39.43 | 39.43 | 39.43 Mbytes
Step TotEng KinEng PotEng E_bond E_angle E_vdwl E_coul E_long Press Volume Density c_TATOM c_TDRUDE[1] c_TDRUDE[2]
0 2707.9246 1420.362 1287.5627 1474.8647 30.734202 -461.31379 540022.21 -539778.93 -9408.4231 42875 0.48077588 298 294.34401 6.9394757
50 1418.4664 1249.502 168.96443 309.99547 480.01554 -495.04808 539716.53 -540017.02 3344.2851 42776.949 0.48187789 252.07241 255.56807 12.841403
100 1391.631 1287.2692 104.36186 333.86636 427.43215 -456.38003 539638.17 -540023.61 -3962.1786 42668.283 0.48310512 265.74273 268.23371 3.3477444
150 1404.7861 1334.1949 70.591169 285.69315 480.49154 -340.69071 539518.01 -540045.75 -2082.9182 42395.713 0.48621109 274.19438 278.9463 1.6008228
200 1425.8199 1378.4982 47.321668 460.35813 459.3409 -415.44782 539411.67 -540044.6 6205.8012 42117.045 0.48942811 284.32248 288.46225 1.1474711
250 1451.5312 1408.0791 43.452124 376.36315 598.2256 -302.5574 539235.03 -540053.8 -3523.5224 41843.302 0.49263001 290.29735 294.90585 0.66499891
300 1460.3317 1532.45 -72.118296 394.63946 561.19451 -334.56109 539192.83 -540050.12 4984.7253 41569.486 0.49587493 316.03909 320.8862 0.85911251
350 1465.3219 1564.1529 -98.830958 418.42916 585.9688 -296.62501 539077.24 -540036.78 740.62275 41386.895 0.49806264 322.55774 327.69261 0.54088355
400 1446.0613 1528.77 -82.708678 414.47172 582.41623 -279.89537 539064.81 -540045.35 -1758.0709 41241.935 0.49981326 314.95112 320.11103 0.86624717
450 1411.6647 1569.3338 -157.66914 380.54196 668.12309 -290.80881 538968.32 -540032.63 3122.7219 41154.153 0.50087936 323.38401 328.46391 1.1709084
500 1366.0173 1653.6398 -287.62246 433.22615 598.36853 -281.95956 538837.44 -540025.77 -2920.9274 41080.823 0.50177344 340.53461 346.03932 1.3737011
Loop time of 21.2526 on 1 procs for 500 steps with 2800 atoms
Performance: 2.033 ns/day, 11.807 hours/ns, 23.526 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.0143 | 5.0143 | 5.0143 | 0.0 | 23.59
Bond | 0.34285 | 0.34285 | 0.34285 | 0.0 | 1.61
Kspace | 7.6454 | 7.6454 | 7.6454 | 0.0 | 35.97
Neigh | 0.36282 | 0.36282 | 0.36282 | 0.0 | 1.71
Comm | 0.035159 | 0.035159 | 0.035159 | 0.0 | 0.17
Output | 0.00068069 | 0.00068069 | 0.00068069 | 0.0 | 0.00
Modify | 7.8451 | 7.8451 | 7.8451 | 0.0 | 36.91
Other | | 0.006337 | | | 0.03
Nlocal: 2800.00 ave 2800 max 2800 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 6816.00 ave 6816 max 6816 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 446966.0 ave 446966 max 446966 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 446966
Ave neighs/atom = 159.63071
Ave special neighs/atom = 11.714286
Neighbor list builds = 28
Dangerous builds = 0
Total wall time: 0:00:21

View File

@ -0,0 +1,200 @@
LAMMPS (24 Aug 2020)
using 1 OpenMP thread(s) per MPI task
# created by fftool
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
special_bonds lj/coul 0.0 0.0 0.5
pair_style hybrid/overlay lj/cut/coul/long 8.0 8.0 thole 2.600 8.0 coul/tt 4 8.0
pair_modify tail yes
kspace_style pppm 1.0e-5
read_data data.ethylene_glycol
Reading data file ...
orthogonal box = (0.0000000 0.0000000 0.0000000) to (35.000000 35.000000 35.000000)
1 by 2 by 2 MPI processor grid
reading atoms ...
2800 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
9 = max dihedrals/atom
reading bonds ...
2600 bonds
reading angles ...
2800 angles
reading dihedrals ...
3000 dihedrals
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0.0 0.0 0.5
special bond factors coul: 0.0 0.0 0.5
5 = max # of 1-2 neighbors
6 = max # of 1-3 neighbors
10 = max # of 1-4 neighbors
13 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.010 seconds
pair_coeff 1 1 lj/cut/coul/long 0.057289 3.500000 # CTO CTO ~
pair_coeff 1 2 lj/cut/coul/long 0.091945 3.304542 # CTO OHG ~
pair_coeff 1 3 lj/cut/coul/long 0.038625 2.958040 # CTO H1O ~
pair_coeff 1 4 lj/cut/coul/long 0.000000 0.000000 # CTO HOG ~
pair_coeff 2 2 lj/cut/coul/long 0.147565 3.120000 # OHG OHG ~
pair_coeff 2 3 lj/cut/coul/long 0.061990 2.792848 # OHG H1O ~
pair_coeff 2 4 lj/cut/coul/long 0.000000 0.000000 # OHG HOG ~
pair_coeff 3 3 lj/cut/coul/long 0.026041 2.500000 # H1O H1O ~
pair_coeff 3 4 lj/cut/coul/long 0.000000 0.000000 # H1O HOG ~
pair_coeff 4 4 lj/cut/coul/long 0.000000 0.000000 # HOG HOG ~
pair_coeff * 5* lj/cut/coul/long 0.000000 0.000000
pair_coeff 1 1 thole 1.662
pair_coeff 1 2 thole 1.561
pair_coeff 1 5 thole 1.662
pair_coeff 1 6 thole 1.561
pair_coeff 2 2 thole 1.467
pair_coeff 2 5 thole 1.561
pair_coeff 2 6 thole 1.467
pair_coeff 5 5 thole 1.662
pair_coeff 5 6 thole 1.561
pair_coeff 6 6 thole 1.467
pair_coeff 2 4 coul/tt 4.5 1.0
pair_coeff 4 6 coul/tt 4.5 1.0
pair_coeff 1 4 coul/tt 4.5 1.0
pair_coeff 4 5 coul/tt 4.5 1.0
group ATOMS type 1 2 3 4
2000 atoms in group ATOMS
group CORES type 1 2
800 atoms in group CORES
group DRUDES type 5 6
800 atoms in group DRUDES
fix DRUDE all drude C C N N D D
fix SHAKE ATOMS shake 0.0001 20 0 b 3 4
400 = # of size 2 clusters
400 = # of size 3 clusters
0 = # of size 4 clusters
0 = # of frozen angles
find clusters CPU = 0.000 seconds
neighbor 2.0 bin
timestep 1.0
variable TK equal 298.0
variable TDK equal 1.0
variable PBAR equal 1.0
comm_modify vel yes
velocity ATOMS create ${TK} 12345
velocity ATOMS create 298 12345
compute TATOM ATOMS temp
compute TDRUDE all temp/drude
fix DTDIR all drude/transform/direct
fix TSTAT ATOMS npt temp ${TK} ${TK} 200 iso ${PBAR} ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 ${TK} 200 iso ${PBAR} ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 298 200 iso ${PBAR} ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 298 200 iso 1 ${PBAR} 1000
fix TSTAT ATOMS npt temp 298 298 200 iso 1 1 1000
fix_modify TSTAT temp TATOM press thermo_press
WARNING: Temperature for fix modify is not for group all (src/fix_nh.cpp:1428)
fix TSTDR DRUDES nvt temp ${TDK} ${TDK} 50
fix TSTDR DRUDES nvt temp 1 ${TDK} 50
fix TSTDR DRUDES nvt temp 1 1 50
fix DTINV all drude/transform/inverse
fix ICECUBE all momentum 1000 linear 1 1 1
thermo_style custom step etotal ke pe ebond eangle evdwl ecoul elong press vol density c_TATOM c_TDRUDE[1] c_TDRUDE[2]
thermo 50
run 500
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:328)
G vector (1/distance) = 0.41206781
grid = 54 54 54
stencil order = 5
estimated absolute RMS force accuracy = 0.0040479865
estimated relative force accuracy = 1.2190391e-05
using double precision FFTW3
3d grid and FFT values/proc = 81648 40824
Rebuild special list taking Drude particles into account
Old max number of 1-2 to 1-4 neighbors: 13
New max number of 1-2 to 1-4 neighbors: 13 (+0)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10
ghost atom cutoff = 10
binsize = 5, bins = 7 7 7
3 neighbor lists, perpetual/occasional/extra = 3 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair thole, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(3) pair coul/tt, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 21.28 | 21.36 | 21.43 Mbytes
Step TotEng KinEng PotEng E_bond E_angle E_vdwl E_coul E_long Press Volume Density c_TATOM c_TDRUDE[1] c_TDRUDE[2]
0 2707.9246 1420.362 1287.5627 1474.8647 30.734202 -461.31379 540022.21 -539778.93 -9408.4231 42875 0.48077588 298 294.34401 6.9394757
50 1418.4665 1249.502 168.96445 309.99547 480.01554 -495.04808 539716.53 -540017.02 3344.2849 42776.949 0.48187789 252.07241 255.56807 12.841403
100 1391.631 1287.2692 104.36189 333.86638 427.43215 -456.38003 539638.17 -540023.61 -3962.1786 42668.283 0.48310512 265.74272 268.23371 3.3477454
150 1404.786 1334.1949 70.59104 285.69317 480.49154 -340.6907 539518.01 -540045.75 -2082.9181 42395.713 0.48621109 274.19439 278.9463 1.6008224
200 1425.8201 1378.4982 47.32185 460.35819 459.34089 -415.44783 539411.67 -540044.6 6205.8014 42117.045 0.48942811 284.32247 288.46225 1.1474763
250 1451.5312 1408.0791 43.45212 376.3631 598.22557 -302.55737 539235.03 -540053.8 -3523.5218 41843.302 0.49263001 290.29734 294.90584 0.66500219
300 1460.3316 1532.4498 -72.118238 394.63949 561.1946 -334.56112 539192.83 -540050.12 4984.7226 41569.486 0.49587493 316.03906 320.88617 0.85911374
350 1465.322 1564.1496 -98.827597 418.43029 585.96947 -296.62447 539077.24 -540036.78 740.63702 41386.895 0.49806264 322.55705 327.69192 0.54088385
400 1446.0613 1528.7739 -82.712576 414.46626 582.41806 -279.89551 539064.81 -540045.35 -1758.0964 41241.935 0.49981326 314.95192 320.11185 0.86625313
450 1411.6627 1569.3001 -157.63739 380.53709 668.13885 -290.80678 538968.36 -540032.62 3122.8773 41154.153 0.50087936 323.37686 328.45689 1.170821
500 1366.0164 1653.6078 -287.59139 433.23854 598.39638 -281.95272 538837.43 -540025.77 -2921.1108 41080.827 0.50177339 340.52802 346.03257 1.3738053
Loop time of 8.54338 on 4 procs for 500 steps with 2800 atoms
Performance: 5.057 ns/day, 4.746 hours/ns, 58.525 timesteps/s
96.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.1432 | 1.2405 | 1.3019 | 5.3 | 14.52
Bond | 0.085337 | 0.091544 | 0.097553 | 1.4 | 1.07
Kspace | 4.3889 | 4.4519 | 4.5415 | 2.6 | 52.11
Neigh | 0.098468 | 0.098565 | 0.098639 | 0.0 | 1.15
Comm | 0.090891 | 0.096061 | 0.102 | 1.3 | 1.12
Output | 0.00051045 | 0.00081128 | 0.0016632 | 0.0 | 0.01
Modify | 2.549 | 2.5525 | 2.5554 | 0.1 | 29.88
Other | | 0.01154 | | | 0.14
Nlocal: 700.000 ave 721 max 665 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Nghost: 4229.75 ave 4342 max 4122 min
Histogram: 1 1 0 0 0 0 0 1 0 1
Neighs: 111742.0 ave 117831 max 101992 min
Histogram: 1 0 0 0 0 1 0 0 0 2
Total # of neighbors = 446967
Ave neighs/atom = 159.63107
Ave special neighs/atom = 11.714286
Neighbor list builds = 28
Dangerous builds = 0
Total wall time: 0:00:08

View File

@ -16,8 +16,8 @@
! and library.h using the ISO_C_BINDING module of the Fortran compiler.
!
! Based on the LAMMPS Fortran 2003 module contributed by:
! Karl D. Hammond <karlh@ugcs.caltech.edu>
! University of Tennessee, Knoxville (USA), 2012
! Karl D. Hammond <hammondkd@missouri.edu>
! University of Missouri, 2012-2020
!
! The Fortran module tries to follow the API of the C-library interface
! closely, but like the Python wrapper it employs an object oriented
@ -30,7 +30,7 @@
MODULE LIBLAMMPS
USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_ptr, c_null_ptr, c_loc, &
c_int, c_char, c_null_char, c_double
c_int, c_char, c_null_char, c_double, c_size_t, c_f_pointer
IMPLICIT NONE
PRIVATE
@ -114,6 +114,12 @@ MODULE LIBLAMMPS
TYPE(c_ptr), VALUE :: str
END SUBROUTINE lammps_commands_string
FUNCTION lammps_malloc(size) BIND(C, name='malloc')
IMPORT :: c_ptr, c_size_t
INTEGER(c_size_t), value :: size
TYPE(c_ptr) :: lammps_malloc
END FUNCTION lammps_malloc
SUBROUTINE lammps_free(ptr) BIND(C, name='lammps_free')
IMPORT :: c_ptr
TYPE(c_ptr), VALUE :: ptr
@ -267,14 +273,14 @@ CONTAINS
CHARACTER (len=*), INTENT(in) :: f_string
CHARACTER (len=1, kind=c_char), POINTER :: c_string(:)
TYPE(c_ptr) :: ptr
INTEGER :: i, n
INTEGER(c_size_t) :: i, n
n = LEN_TRIM(f_string)
ALLOCATE(c_string(n+1))
ptr = lammps_malloc(n+1)
CALL C_F_POINTER(ptr,c_string,[1])
DO i=1,n
c_string(i) = f_string(i:i)
END DO
c_string(n+1) = c_null_char
ptr = c_loc(c_string(1))
END FUNCTION f2c_string
END MODULE LIBLAMMPS

View File

@ -7,6 +7,7 @@ features:
using Langevin or Nosé-Hoover thermostats
* computation of the atom and dipole temperatures
* damping induced dipole interactions using Thole's function
* charge-dipole damping using Tang-Toennies function
See the file doc/drude_tutorial.html for getting started.
@ -18,4 +19,5 @@ The person who created this package is Alain Dequidt at the
Chemistry Institute of Clermont-Ferrand, Clermont University, France
(alain.dequidt at uca.fr). Contact him directly if you have questions.
Co-authors: Julien Devémy, Agilio Padua.
Contributors: Kateryna Goloviznina, Zheng Gong.

View File

@ -0,0 +1,446 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "pair_coul_tt.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "fix.h"
#include "fix_drude.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "modify.h"
#include "error.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairCoulTT::PairCoulTT(LAMMPS *lmp) : Pair(lmp) {
fix_drude = nullptr;
}
/* ---------------------------------------------------------------------- */
PairCoulTT::~PairCoulTT()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(b);
memory->destroy(c);
memory->destroy(ntt);
memory->destroy(cut);
memory->destroy(scale);
}
}
/* ---------------------------------------------------------------------- */
void PairCoulTT::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qi,qj,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double r,rsq,r2inv,rinv,factor_coul;
int *ilist,*jlist,*numneigh,**firstneigh;
double factor_f,factor_e;
int di,dj;
double dcoul;
double beta, gamma, betaprime, gammaprime, gammatmp;
ecoul = 0.0;
ev_init(eflag,vflag);
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
int *drudetype = fix_drude->drudetype;
tagint *drudeid = fix_drude->drudeid;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qi = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
if (drudetype[type[i]] == drudetype[type[j]] && drudetype[type[j]] != CORE_TYPE)
continue;
qj = q[j];
if (drudetype[type[i]] == CORE_TYPE) {
di = domain->closest_image(i, atom->map(drudeid[i]));
if (di == j)
continue;
switch (drudetype[type[j]]) {
case DRUDE_TYPE:
qi = q[i]+q[di];
break;
case NOPOL_TYPE:
qi = -q[di];
break;
}
}
if (drudetype[type[j]] == CORE_TYPE) {
dj = domain->closest_image(j, atom->map(drudeid[j]));
if (dj == i)
continue;
switch (drudetype[type[i]]) {
case DRUDE_TYPE:
qj = q[j]+q[dj];
break;
case NOPOL_TYPE:
qj = -q[dj];
break;
}
}
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r = sqrt(rsq);
beta = c[itype][jtype] * exp(-b[itype][jtype] * r);
betaprime = -b[itype][jtype] * beta;
gamma = 1.0 + b[itype][jtype] * r;
gammaprime = b[itype][jtype];
gammatmp = 1.0;
for (int k = 2; k <= ntt[itype][jtype]; k++) {
gammatmp *= b[itype][jtype] * r / k;
gamma += gammatmp * b[itype][jtype] * r;
gammaprime += gammatmp * b[itype][jtype] * k;
}
if (drudetype[type[i]] == CORE_TYPE && drudetype[type[j]] == CORE_TYPE)
dcoul = qqrd2e * ( -(q[i]+q[di])*q[dj] - q[di]*(q[j]+q[dj]) ) * scale[itype][jtype] * rinv;
else
dcoul = qqrd2e * qi * qj *scale[itype][jtype] * rinv;
factor_f = (-beta*gamma + r*betaprime*gamma + r*beta*gammaprime)*factor_coul;
if(eflag) factor_e = - beta*gamma*factor_coul;
fpair = factor_f * dcoul * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag)
ecoul = factor_e * dcoul;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairCoulTT::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(scale,n+1,n+1,"pair:scale");
memory->create(b, n + 1, n + 1, "pair:b");
memory->create(c, n + 1, n + 1, "pair:c");
memory->create(ntt, n + 1, n + 1, "pair:ntt");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairCoulTT::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR, "Illegal pair_style command");
n_global = utils::inumeric(FLERR,arg[0],false,lmp);
cut_global = utils::numeric(FLERR,arg[1],false,lmp);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) {
ntt[i][j] = n_global;
cut[i][j] = cut_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairCoulTT::coeff(int narg, char **arg)
{
if (narg < 3 || narg > 6)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
double b_one = utils::numeric(FLERR,arg[2],false,lmp);
double c_one = utils::numeric(FLERR,arg[3],false,lmp);
int n_one = n_global;
double cut_one = cut_global;
if (narg >= 5) n_one = utils::inumeric(FLERR,arg[4],false,lmp);
if (narg == 6) cut_one = utils::numeric(FLERR,arg[5],false,lmp);
if (n_one > n_global)
error->all(FLERR, "Incorrect coefficients for pair style coul/tt: n should not be larger than global setting");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
b[i][j] = b_one;
c[i][j] = c_one;
ntt[i][j] = n_one;
cut[i][j] = cut_one;
scale[i][j] = 1.0;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairCoulTT::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style coul/tt requires atom attribute q");
int ifix;
for (ifix = 0; ifix < modify->nfix; ifix++)
if (utils::strmatch(modify->fix[ifix]->style,"^drude")) break;
if (ifix == modify->nfix) error->all(FLERR, "Pair coul/tt requires fix drude");
fix_drude = (FixDrude *) modify->fix[ifix];
neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairCoulTT::init_one(int i, int j)
{
if (setflag[i][j] == 0)
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
b[j][i] = b[i][j];
c[j][i] = c[i][j];
ntt[j][i] = ntt[i][j];
scale[j][i] = scale[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulTT::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&b[i][j], sizeof(double), 1, fp);
fwrite(&c[i][j], sizeof(double), 1, fp);
fwrite(&ntt[i][j], sizeof(int), 1, fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulTT::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR,&b[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&c[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&ntt[i][j],sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
}
MPI_Bcast(&b[i][j], 1, MPI_DOUBLE,0,world);
MPI_Bcast(&c[i][j], 1, MPI_DOUBLE,0,world);
MPI_Bcast(&ntt[i][j], 1, MPI_INT,0,world);
MPI_Bcast(&cut[i][j],1, MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairCoulTT::write_restart_settings(FILE *fp)
{
fwrite(&n_global, sizeof(int), 1, fp);
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairCoulTT::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR,&n_global,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
}
MPI_Bcast(&n_global, 1, MPI_INT, 0, world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairCoulTT::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double /*factor_lj*/,
double &fforce)
{
double r2inv,rinv,r,phicoul;
double qi,qj,factor_f,factor_e,dcoul;
double beta, betaprime, gamma, gammaprime, gammatmp;
// single() is used in Pair::write_file where there is no information about topology
// it corresponds to pair_write command in input file
// charges qi and qj are defined by the user (or 1.0 by defaut)
qi = atom->q[i];
qj = atom->q[j];
r2inv = 1.0/rsq;
fforce = phicoul = 0.0;
if (rsq < cutsq[itype][jtype]) {
rinv = sqrt(r2inv);
r = sqrt(rsq);
beta = c[itype][jtype] * exp(-b[itype][jtype] * r);
betaprime = -b[itype][jtype] * beta;
gamma = 1 + b[itype][jtype] * r;
gammaprime = b[itype][jtype];
gammatmp = 1;
for (int k = 2; k <= ntt[itype][jtype]; k++) {
gammatmp *= b[itype][jtype] * r / k;
gamma += gammatmp * b[itype][jtype] * r;
gammaprime += gammatmp * b[itype][jtype] * k;
}
dcoul = force->qqrd2e * qi * qj * scale[itype][jtype] * rinv;
factor_f = (-beta*gamma+r*betaprime*gamma+r*beta*gammaprime)*factor_coul;
fforce = factor_f * dcoul * r2inv;
factor_e = - beta*gamma*factor_coul;
phicoul = factor_e * dcoul;
}
return phicoul;
}
/* ---------------------------------------------------------------------- */
void *PairCoulTT::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"scale") == 0) return (void *) scale;
if (strcmp(str,"b") == 0) return (void *) b;
if (strcmp(str,"c") == 0) return (void *) c;
if (strcmp(str,"ntt") == 0) return (void *) ntt;
return nullptr;
}

View File

@ -0,0 +1,75 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(coul/tt,PairCoulTT)
#else
#ifndef LMP_PAIR_COULTT_H
#define LMP_PAIR_COULTT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairCoulTT: public Pair {
public:
PairCoulTT(class LAMMPS *);
virtual ~PairCoulTT();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
void *extract(const char *, int &);
protected:
int n_global;
double cut_global;
double **cut,**scale;
double **b,**c;
int **ntt;
class FixDrude * fix_drude;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style thole requires atom attribute q
The atom style defined does not have these attributes.
*/