lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3539 from akohlmey/collected-small-changes 

Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer 2022-11-30 13:45:05 -05:00 committed by GitHub
parent 63b2d2eec7 72140e2608
commit e76e864152
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
36 changed files with 813 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.github/CODEOWNERS vendored
View File

@ -37,6 +37,7 @@ src/MESONT/* @iafoss
src/ML-HDNNP/* @singraber
src/ML-IAP/* @athomps
src/ML-PACE/* @yury-lysogorskiy
src/ML-POD/* @exapde @rohskopf
src/MOFFF/* @hheenen
src/MOLFILE/* @akohlmey
src/NETCDF/* @pastewka

8
doc/src/Developer_code_design.rst
View File

@ -78,7 +78,7 @@ LAMMPS makes extensive use of the object oriented programming (OOP)
principles of *compositing* and *inheritance*. Classes like the
``LAMMPS`` class are a **composite** containing pointers to instances
of other classes like ``Atom``, ``Comm``, ``Force``, ``Neighbor``,
``Modify``, and so on. Each of these classes implement certain
``Modify``, and so on. Each of these classes implements certain
functionality by storing and manipulating data related to the
simulation and providing member functions that trigger certain
actions. Some of those classes like ``Force`` are themselves
@ -87,9 +87,9 @@ interactions. Similarly the ``Modify`` class contains a list of
``Fix`` and ``Compute`` classes. If the input commands that
correspond to these classes include the word *style*, then LAMMPS
stores only a single instance of that class. E.g. *atom_style*,
*comm_style*, *pair_style*, *bond_style*. It the input command does
not include the word *style*, there can be many instances of that
class defined. E.g. *region*, *fix*, *compute*, *dump*.
*comm_style*, *pair_style*, *bond_style*. If the input command does
**not** include the word *style*, then there may be many instances of
that class defined, for example *region*, *fix*, *compute*, *dump*.
**Inheritance** enables creation of *derived* classes that can share
common functionality in their base class while providing a consistent

2
doc/src/compute_pair_local.rst
View File

@ -59,7 +59,7 @@ commands.
The value *dist* is the distance between the pair of atoms.
The values *dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the
*distance* between the pair of atoms. This value is always the
distance from the atom of lower to the one with the higher id.
distance from the atom of higher to the one with the lower atom ID.
The value *eng* is the interaction energy for the pair of atoms.

6
doc/src/fix_ipi.rst
View File

@ -90,6 +90,12 @@ coordinates are transferred. However, one could use this strategy to
define an external potential acting on the atoms that are moved by
i-PI.
Since the i-PI code uses atomic units internally, this fix needs to
convert LAMMPS data to and from its :doc:`specified units <units>`
accordingly when communicating with i-PI. This is not possible for
reduced units ("units lj") and thus *fix ipi* will stop with an error in
this case.
This fix is part of the MISC package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.

2
doc/src/fix_pimd.rst
View File

@ -156,6 +156,8 @@ This fix is part of the REPLICA package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
Fix pid cannot be used with :doc:`lj units <units>`.
A PIMD simulation can be initialized with a single data file read via
the :doc:`read_data <read_data>` command. However, this means all
quasi-beads in a ring polymer will have identical positions and

10
doc/src/pair_sw.rst
View File

@ -176,9 +176,13 @@ are placeholders for atom types that will be used with other potentials.
.. note::
When the *threebody off* keyword is used, multiple pair_coeff commands may
be used to specific the pairs of atoms which don't require three-body term.
In these cases, the first 2 arguments are not required to be \* \*.
When the *threebody off* keyword is used, multiple pair_coeff
commands may be used to specific the pairs of atoms which don't
require three-body term. In these cases, the first 2 arguments are
not required to be \* \*, the potential parameter file is only read
by the first :doc:`pair_coeff command <pair_coeff>` and the element
to atom type mappings must be consistent across all *pair_coeff*
statements. If not LAMMPS will abort with an error.
Stillinger-Weber files in the *potentials* directory of the LAMMPS
distribution have a ".sw" suffix. Lines that are not blank or

16
lib/linalg/README
View File

@ -1,13 +1,16 @@
This directory has BLAS and LAPACK files needed by the ATC and
AWPMD packages, and possibly by other packages in the future.
This directory has generic BLAS and LAPACK source files needed by the
ATC, AWPMD, ELECTRODE, LATTE, and ML-POD packages (and possibly by other
packages) in the future that can be used instead of platform or vendor
optimized BLAS/LAPACK library.
Note that this is an *incomplete* subset of full BLAS/LAPACK.
The files correspond to LAPACK version 3.11.0.
You should only need to build and use the library in this directory if
you want to build LAMMPS with the ATC and/or AWPMD packages
AND you do not have any other suitable BLAS and LAPACK libraries
installed on your system. E.g. ATLAS, GOTO-BLAS, OpenBLAS, ACML, or
MKL.
you want to build LAMMPS with one of the listed packages AND you do not
have any other suitable BLAS and LAPACK libraries installed on your
system (like ATLAS, GOTO-BLAS, OpenBLAS, ACML, MKL and so on).
You can type "make lib-linalg" from the src directory to see help on
how to build this library via make commands, or you can do the same
@ -27,3 +30,4 @@ liblinalg.a the library LAMMPS will link against
You can then include this library and its path in the Makefile.lammps
file of any packages that need it. As an example, see the
lib/atc/Makefile.lammps.linalg file.

26
lib/linalg/dgelss.f
View File

@ -254,11 +254,11 @@
*
* Compute space needed for DGEQRF
CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
LWORK_DGEQRF=DUM(1)
LWORK_DGEQRF = INT( DUM(1) )
* Compute space needed for DORMQR
CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B,
$ LDB, DUM(1), -1, INFO )
LWORK_DORMQR=DUM(1)
LWORK_DORMQR = INT( DUM(1) )
MM = N
MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF )
MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR )
@ -273,15 +273,15 @@
* Compute space needed for DGEBRD
CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, INFO )
LWORK_DGEBRD=DUM(1)
LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORMBR
CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1),
$ B, LDB, DUM(1), -1, INFO )
LWORK_DORMBR=DUM(1)
LWORK_DORMBR = INT( DUM(1) )
* Compute space needed for DORGBR
CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
$ DUM(1), -1, INFO )
LWORK_DORGBR=DUM(1)
LWORK_DORGBR = INT( DUM(1) )
* Compute total workspace needed
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR )
@ -305,23 +305,23 @@
* Compute space needed for DGELQF
CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1),
$ -1, INFO )
LWORK_DGELQF=DUM(1)
LWORK_DGELQF = INT( DUM(1) )
* Compute space needed for DGEBRD
CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, INFO )
LWORK_DGEBRD=DUM(1)
LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORMBR
CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA,
$ DUM(1), B, LDB, DUM(1), -1, INFO )
LWORK_DORMBR=DUM(1)
LWORK_DORMBR = INT( DUM(1) )
* Compute space needed for DORGBR
CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1),
$ DUM(1), -1, INFO )
LWORK_DORGBR=DUM(1)
LWORK_DORGBR = INT( DUM(1) )
* Compute space needed for DORMLQ
CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1),
$ B, LDB, DUM(1), -1, INFO )
LWORK_DORMLQ=DUM(1)
LWORK_DORMLQ = INT( DUM(1) )
* Compute total workspace needed
MAXWRK = M + LWORK_DGELQF
MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD )
@ -341,15 +341,15 @@
* Compute space needed for DGEBRD
CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
$ DUM(1), DUM(1), -1, INFO )
LWORK_DGEBRD=DUM(1)
LWORK_DGEBRD = INT( DUM(1) )
* Compute space needed for DORMBR
CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA,
$ DUM(1), B, LDB, DUM(1), -1, INFO )
LWORK_DORMBR=DUM(1)
LWORK_DORMBR = INT( DUM(1) )
* Compute space needed for DORGBR
CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1),
$ DUM(1), -1, INFO )
LWORK_DORGBR=DUM(1)
LWORK_DORGBR = INT( DUM(1) )
MAXWRK = 3*M + LWORK_DGEBRD
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR )
MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR )

7
lib/linalg/dlaed4.f
View File

@ -328,9 +328,12 @@
IF( C.LT.ZERO )
$ C = ABS( C )
IF( C.EQ.ZERO ) THEN
* ETA = B/A
* ETA = B/A
* ETA = RHO - TAU
ETA = DLTUB - TAU
* ETA = DLTUB - TAU
*
* Update proposed by Li, Ren-Cang:
ETA = -W / ( DPSI+DPHI )
ELSE IF( A.GE.ZERO ) THEN
ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE

2
lib/linalg/dlascl.f
View File

@ -272,6 +272,8 @@
ELSE
MUL = CTOC / CFROMC
DONE = .TRUE.
IF (MUL .EQ. ONE)
$ RETURN
END IF
END IF
*

69
lib/linalg/dlatrs.f
View File

@ -264,8 +264,8 @@
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DASUM, DDOT, DLAMCH
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA
@ -304,6 +304,7 @@
*
* Quick return if possible
*
SCALE = ONE
IF( N.EQ.0 )
$ RETURN
*
@ -311,7 +312,6 @@
*
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
@ -343,8 +343,67 @@
IF( TMAX.LE.BIGNUM ) THEN
TSCAL = ONE
ELSE
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
*
* Avoid NaN generation if entries in CNORM exceed the
* overflow threshold
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
* Case 1: All entries in CNORM are valid floating-point numbers
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
ELSE
* Case 2: At least one column norm of A cannot be represented
* as floating-point number. Find the offdiagonal entry A( I, J )
* with the largest absolute value. If this entry is not +/- Infinity,
* use this value as TSCAL.
TMAX = ZERO
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO J = 2, N
TMAX = MAX( DLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
$ TMAX )
END DO
ELSE
*
* A is lower triangular.
*
DO J = 1, N - 1
TMAX = MAX( DLANGE( 'M', N-J, 1, A( J+1, J ), 1,
$ SUMJ ), TMAX )
END DO
END IF
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
TSCAL = ONE / ( SMLNUM*TMAX )
DO J = 1, N
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
CNORM( J ) = CNORM( J )*TSCAL
ELSE
* Recompute the 1-norm without introducing Infinity
* in the summation
CNORM( J ) = ZERO
IF( UPPER ) THEN
DO I = 1, J - 1
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
ELSE
DO I = J + 1, N
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
END IF
END IF
END DO
ELSE
* At least one entry of A is not a valid floating-point entry.
* Rely on TRSV to propagate Inf and NaN.
CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
RETURN
END IF
END IF
END IF
*
* Compute a bound on the computed solution vector to see if the

2
lib/linalg/dorgbr.f
View File

@ -232,7 +232,7 @@
END IF
END IF
END IF
LWKOPT = WORK( 1 )
LWKOPT = INT( WORK( 1 ) )
LWKOPT = MAX (LWKOPT, MN)
END IF
*

190
lib/linalg/dposv.f Normal file
View File

@ -0,0 +1,190 @@
*> \brief <b> DPOSV computes the solution to system of linear equations A * X = B for PO matrices</b>
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DPOSV + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dposv.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dposv.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dposv.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DPOSV computes the solution to a real system of linear equations
*> A * X = B,
*> where A is an N-by-N symmetric positive definite matrix and X and B
*> are N-by-NRHS matrices.
*>
*> The Cholesky decomposition is used to factor A as
*> A = U**T* U, if UPLO = 'U', or
*> A = L * L**T, if UPLO = 'L',
*> where U is an upper triangular matrix and L is a lower triangular
*> matrix. The factored form of A is then used to solve the system of
*> equations A * X = B.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> = 'U': Upper triangle of A is stored;
*> = 'L': Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of linear equations, i.e., the order of the
*> matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*> NRHS is INTEGER
*> The number of right hand sides, i.e., the number of columns
*> of the matrix B. NRHS >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
*> N-by-N upper triangular part of A contains the upper
*> triangular part of the matrix A, and the strictly lower
*> triangular part of A is not referenced. If UPLO = 'L', the
*> leading N-by-N lower triangular part of A contains the lower
*> triangular part of the matrix A, and the strictly upper
*> triangular part of A is not referenced.
*>
*> On exit, if INFO = 0, the factor U or L from the Cholesky
*> factorization A = U**T*U or A = L*L**T.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
*> On entry, the N-by-NRHS right hand side matrix B.
*> On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> > 0: if INFO = i, the leading minor of order i of A is not
*> positive definite, so the factorization could not be
*> completed, and the solution has not been computed.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup doublePOsolve
*
* =====================================================================
SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DPOTRF, DPOTRS, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( NRHS.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -7
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DPOSV ', -INFO )
RETURN
END IF
*
* Compute the Cholesky factorization A = U**T*U or A = L*L**T.
*
CALL DPOTRF( UPLO, N, A, LDA, INFO )
IF( INFO.EQ.0 ) THEN
*
* Solve the system A*X = B, overwriting B with X.
*
CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
END IF
RETURN
*
* End of DPOSV
*
END

201
lib/linalg/dpotrs.f Normal file
View File

@ -0,0 +1,201 @@
*> \brief \b DPOTRS
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DPOTRS + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpotrs.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpotrs.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpotrs.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DPOTRS solves a system of linear equations A*X = B with a symmetric
*> positive definite matrix A using the Cholesky factorization
*> A = U**T*U or A = L*L**T computed by DPOTRF.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> = 'U': Upper triangle of A is stored;
*> = 'L': Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The order of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*> NRHS is INTEGER
*> The number of right hand sides, i.e., the number of columns
*> of the matrix B. NRHS >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> The triangular factor U or L from the Cholesky factorization
*> A = U**T*U or A = L*L**T, as computed by DPOTRF.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
*> On entry, the right hand side matrix B.
*> On exit, the solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup doublePOcomputational
*
* =====================================================================
SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
*
* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL UPPER
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DTRSM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
UPPER = LSAME( UPLO, 'U' )
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( NRHS.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -7
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DPOTRS', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 .OR. NRHS.EQ.0 )
$ RETURN
*
IF( UPPER ) THEN
*
* Solve A*X = B where A = U**T *U.
*
* Solve U**T *X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
$ ONE, A, LDA, B, LDB )
*
* Solve U*X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
$ NRHS, ONE, A, LDA, B, LDB )
ELSE
*
* Solve A*X = B where A = L*L**T.
*
* Solve L*X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N,
$ NRHS, ONE, A, LDA, B, LDB )
*
* Solve L**T *X = B, overwriting B with X.
*
CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', N, NRHS,
$ ONE, A, LDA, B, LDB )
END IF
*
RETURN
*
* End of DPOTRS
*
END

5
lib/linalg/dscal.f
View File

@ -93,11 +93,14 @@
*
* .. Local Scalars ..
INTEGER I,M,MP1,NINCX
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER (ONE=1.0D+0)
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1

2
lib/linalg/dsyevd.f
View File

@ -257,7 +257,7 @@
LWMIN = 2*N + 1
END IF
LOPT = MAX( LWMIN, 2*N +
$ ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
$ N*ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
LIOPT = LIWMIN
END IF
WORK( 1 ) = LOPT

4
lib/linalg/dsygvd.f
View File

@ -330,8 +330,8 @@
CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
$ INFO )
LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
LOPT = INT( MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) )
LIOPT = INT( MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) )
*
IF( WANTZ .AND. INFO.EQ.0 ) THEN
*

2
lib/linalg/ieeeck.f
View File

@ -41,7 +41,7 @@
*> \param[in] ISPEC
*> \verbatim
*> ISPEC is INTEGER
*> Specifies whether to test just for inifinity arithmetic
*> Specifies whether to test just for infinity arithmetic
*> or whether to test for infinity and NaN arithmetic.
*> = 0: Verify infinity arithmetic only.
*> = 1: Verify infinity and NaN arithmetic.

15
lib/linalg/ilaenv.f
View File

@ -469,6 +469,15 @@
ELSE
NB = 64
END IF
ELSE IF( C3.EQ.'SYL' ) THEN
* The upper bound is to prevent overly aggressive scaling.
IF( SNAME ) THEN
NB = MIN( MAX( 48, INT( ( MIN( N1, N2 ) * 16 ) / 100) ),
$ 240 )
ELSE
NB = MIN( MAX( 24, INT( ( MIN( N1, N2 ) * 8 ) / 100) ),
$ 80 )
END IF
END IF
ELSE IF( C2.EQ.'LA' ) THEN
IF( C3.EQ.'UUM' ) THEN
@ -477,6 +486,12 @@
ELSE
NB = 64
END IF
ELSE IF( C3.EQ.'TRS' ) THEN
IF( SNAME ) THEN
NB = 32
ELSE
NB = 32
END IF
END IF
ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
IF( C3.EQ.'EBZ' ) THEN

11
lib/linalg/zdscal.f
View File

@ -92,17 +92,20 @@
*
* .. Local Scalars ..
INTEGER I,NINCX
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER (ONE=1.0D+0)
* ..
* .. Intrinsic Functions ..
INTRINSIC DCMPLX
INTRINSIC DBLE, DCMPLX, DIMAG
* ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1
*
DO I = 1,N
ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
ZX(I) = DCMPLX(DA*DBLE(ZX(I)),DA*DIMAG(ZX(I)))
END DO
ELSE
*
@ -110,7 +113,7 @@
*
NINCX = N*INCX
DO I = 1,NINCX,INCX
ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
ZX(I) = DCMPLX(DA*DBLE(ZX(I)),DA*DIMAG(ZX(I)))
END DO
END IF
RETURN

2
lib/linalg/zheevd.f
View File

@ -284,7 +284,7 @@
LIWMIN = 1
END IF
LOPT = MAX( LWMIN, N +
$ ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
$ N*ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
LROPT = LRWMIN
LIOPT = LIWMIN
END IF

2
lib/linalg/zlascl.f
View File

@ -272,6 +272,8 @@
ELSE
MUL = CTOC / CFROMC
DONE = .TRUE.
IF (MUL .EQ. ONE)
$ RETURN
END IF
END IF
*

6
lib/linalg/zscal.f
View File

@ -93,7 +93,11 @@
* .. Local Scalars ..
INTEGER I,NINCX
* ..
IF (N.LE.0 .OR. INCX.LE.0) RETURN
* .. Parameters ..
COMPLEX*16 ONE
PARAMETER (ONE= (1.0D+0,0.0D+0))
* ..
IF (N.LE.0 .OR. INCX.LE.0 .OR. ZA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1

2
src/CG-SPICA/angle_spica.cpp
View File

@ -297,7 +297,7 @@ void AngleSPICA::init_style()
repflag = 0;
for (int i = 1; i <= atom->nangletypes; i++)
if (repscale[i] > 0.0) repflag = 1;
if (repscale && (repscale[i] > 0.0)) repflag = 1;
// set up pointers to access SPICA LJ parameters for 1-3 interactions

56
src/MANYBODY/pair_sw.cpp
View File

@ -46,6 +46,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
centroidstressflag = CENTROID_NOTAVAIL;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
skip_threebody_flag = false;
params_mapped = 0;
params = nullptr;
@ -225,12 +226,12 @@ void PairSW::compute(int eflag, int vflag)
void PairSW::allocate()
{
allocated = 1;
int n = atom->ntypes;
int np1 = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(neighshort,maxshort,"pair:neighshort");
map = new int[n+1];
memory->create(setflag, np1, np1, "pair:setflag");
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(neighshort, maxshort, "pair:neighshort");
map = new int[np1];
}
/* ----------------------------------------------------------------------
@ -262,12 +263,49 @@ void PairSW::coeff(int narg, char **arg)
{
if (!allocated) allocate();
map_element2type(narg-3,arg+3);
// read potential file and set up element maps only once
if (one_coeff || !params_mapped) {
// make certain that the setflag array is always fully initialized
// the sw/intel pair style depends on it
if (!one_coeff) {
for (int i = 0; i <= atom->ntypes; i++) {
for (int j = 0; j <= atom->ntypes; j++) {
setflag[i][j] = 0;
}
}
}
// read potential file and initialize potential parameters
map_element2type(narg-3, arg+3, (one_coeff != 0));
read_file(arg[2]);
setup_params();
// read potential file and initialize potential parameters
read_file(arg[2]);
setup_params();
params_mapped = 1;
}
if (!one_coeff) {
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);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
if (((map[i] >= 0) && (strcmp(arg[i+2], elements[map[i]]) != 0)) ||
((map[i] < 0) && (strcmp(arg[i+2], "NULL") != 0)))
error->all(FLERR, "Must use consistent type to element mappings with threebody off");
if (map[i] < 0) error->all(FLERR, "Must not set pair_coeff mapped to NULL element");
for (int j = MAX(jlo, i); j <= jhi; j++) {
if (((map[j] >= 0) && (strcmp(arg[j+2], elements[map[j]]) != 0)) ||
((map[j] < 0) && (strcmp(arg[j+2], "NULL") != 0)))
error->all(FLERR, "Must use consistent type to element mappings with threebody off");
if (map[j] < 0) error->all(FLERR, "Must not set pair_coeff mapped to NULL element");
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
}
/* ----------------------------------------------------------------------

1
src/MANYBODY/pair_sw.h
View File

@ -54,6 +54,7 @@ class PairSW : public Pair {
int maxshort; // size of short neighbor list array
int *neighshort; // short neighbor list array
int skip_threebody_flag; // whether to run threebody loop
int params_mapped; // whether parameters have been read and mapped to elements
void settings(int, char **) override;
virtual void allocate();

115
src/MISC/fix_ipi.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -47,12 +46,12 @@ using namespace FixConst;
// socket interface
#ifndef _WIN32
#include <unistd.h>
#include <sys/types.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <sys/un.h>
#include <netdb.h>
#include <netinet/in.h>
#include <sys/socket.h>
#include <sys/types.h>
#include <sys/un.h>
#include <unistd.h>
#else
#ifndef WIN32_LEAN_AND_MEAN
#define WIN32_LEAN_AND_MEAN
@ -65,8 +64,7 @@ using namespace FixConst;
/* Utility functions to simplify the interface with POSIX sockets */
static void open_socket(int &sockfd, int inet, int port, char* host,
Error *error)
static void open_socket(int &sockfd, int inet, int port, char *host, Error *error)
/* Opens a socket.
Args:
@ -83,9 +81,9 @@ static void open_socket(int &sockfd, int inet, int port, char* host,
int ai_err;
#ifdef _WIN32
error->one(FLERR,"i-PI socket implementation requires UNIX environment");
error->one(FLERR, "i-PI socket implementation requires UNIX environment");
#else
if (inet>0) { // creates an internet socket
if (inet > 0) { // creates an internet socket
// fetches information on the host
struct addrinfo hints, *res;
@ -96,40 +94,39 @@ static void open_socket(int &sockfd, int inet, int port, char* host,
hints.ai_flags = AI_PASSIVE;
ai_err = getaddrinfo(host, std::to_string(port).c_str(), &hints, &res);
if (ai_err!=0)
error->one(FLERR,"Error fetching host data. Wrong host name?");
if (ai_err != 0) error->one(FLERR, "Error fetching host data. Wrong host name?");
// creates socket
sockfd = socket(res->ai_family, res->ai_socktype, res->ai_protocol);
if (sockfd < 0)
error->one(FLERR,"Error opening socket");
if (sockfd < 0) error->one(FLERR, "Error opening socket");
// makes connection
if (connect(sockfd, res->ai_addr, res->ai_addrlen) < 0)
error->one(FLERR,"Error opening INET socket: wrong port or server unreachable");
error->one(FLERR, "Error opening INET socket: wrong port or server unreachable");
freeaddrinfo(res);
} else { // creates a unix socket
} else { // creates a unix socket
struct sockaddr_un serv_addr;
// fills up details of the socket address
memset(&serv_addr, 0, sizeof(serv_addr));
serv_addr.sun_family = AF_UNIX;
strcpy(serv_addr.sun_path, "/tmp/ipi_");
strcpy(serv_addr.sun_path+9, host);
strcpy(serv_addr.sun_path + 9, host);
// creates the socket
sockfd = socket(AF_UNIX, SOCK_STREAM, 0);
// connects
if (connect(sockfd, (struct sockaddr *) &serv_addr, sizeof(serv_addr)) < 0)
error->one(FLERR,"Error opening UNIX socket: server may not be running "
error->one(FLERR,
"Error opening UNIX socket: server may not be running "
"or the path to the socket unavailable");
}
#endif
}
static void writebuffer(int sockfd, const char *data, int len, Error* error)
static void writebuffer(int sockfd, const char *data, int len, Error *error)
/* Writes to a socket.
Args:
@ -140,14 +137,11 @@ static void writebuffer(int sockfd, const char *data, int len, Error* error)
{
int n;
n = write(sockfd,data,len);
if (n < 0)
error->one(FLERR,"Error writing to socket: broken connection");
n = write(sockfd, data, len);
if (n < 0) error->one(FLERR, "Error writing to socket: broken connection");
}
static void readbuffer(int sockfd, char *data, int len, Error* error)
static void readbuffer(int sockfd, char *data, int len, Error *error)
/* Reads from a socket.
Args:
@ -158,44 +152,52 @@ static void readbuffer(int sockfd, char *data, int len, Error* error)
{
int n, nr;
n = nr = read(sockfd,data,len);
n = nr = read(sockfd, data, len);
while (nr>0 && n<len) {
nr=read(sockfd,&data[n],len-n);
n+=nr;
while (nr > 0 && n < len) {
nr = read(sockfd, &data[n], len - n);
n += nr;
}
if (n == 0)
error->one(FLERR,"Error reading from socket: broken connection");
if (n == 0) error->one(FLERR, "Error reading from socket: broken connection");
}
/* ---------------------------------------------------------------------- */
FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), irregular(nullptr)
FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), irregular(nullptr)
{
/* format for fix:
* fix num group_id ipi host port [unix]
*/
if (strcmp(style,"ipi") != 0 && narg < 5)
error->all(FLERR,"Illegal fix ipi command");
if (narg < 5) utils::missing_cmd_args(FLERR, "fix ipi", error);
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix ipi without atom IDs");
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use fix ipi without atom IDs");
if (atom->tag_consecutive() == 0) error->all(FLERR, "Fix ipi requires consecutive atom IDs");
if (strcmp(update->unit_style, "lj") == 0) error->all(FLERR, "Fix ipi does not support lj units");
if (atom->tag_consecutive() == 0)
error->all(FLERR,"Fix ipi requires consecutive atom IDs");
if (strcmp(arg[1],"all") != 0)
error->warning(FLERR,"Fix ipi always uses group all");
if ((strcmp(arg[1], "all") != 0) && (comm->me == 0))
error->warning(FLERR, "Not using group 'all' with fix ipi can result in undefined behavior");
host = strdup(arg[3]);
port = utils::inumeric(FLERR,arg[4],false,lmp);
port = utils::inumeric(FLERR, arg[4], false, lmp);
inet = ((narg > 5) && (strcmp(arg[5],"unix") == 0) ) ? 0 : 1;
master = (comm->me==0) ? 1 : 0;
// check if forces should be reinitialized and set flag
reset_flag = ((narg > 6 && (strcmp(arg[5],"reset") == 0 )) || ((narg > 5) && (strcmp(arg[5],"reset") == 0)) ) ? 1 : 0;
master = (comm->me == 0) ? 1 : 0;
inet = 1;
reset_flag = 0;
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg], "unix") == 0) {
inet = 0;
++iarg;
} else if (strcmp(arg[iarg], "reset") == 0) {
reset_flag = 1;
++iarg;
} else {
error->all(FLERR, "Unknown fix ipi keyword: {}", arg[iarg]);
}
}
// sanity check
if (inet && ((port <= 1024) || (port > 65536)))
error->all(FLERR, "Invalid port for fix ipi: {}", port);
hasdata = bsize = 0;
@ -223,7 +225,6 @@ FixIPI::~FixIPI()
delete irregular;
}
/* ---------------------------------------------------------------------- */
int FixIPI::setmask()
@ -240,9 +241,11 @@ void FixIPI::init()
{
//only opens socket on master process
if (master) {
if (!socketflag) open_socket(ipisock, inet, port, host, error);
} else ipisock=0;
//! should check for success in socket opening -- but the current open_socket routine dies brutally if unsuccessful
if (!socketflag) open_socket(ipisock, inet, port, host, error);
} else
ipisock = 0;
// TODO: should check for success in socket opening,
// but the current open_socket routine dies brutally if unsuccessful
// tell lammps we have assigned a socket
socketflag = 1;
@ -252,11 +255,13 @@ void FixIPI::init()
kspace_flag = (force->kspace) ? 1 : 0;
// makes sure that neighbor lists are re-built at each step (cannot make assumptions when cycling over beads!)
// makes sure that neighbor lists are re-built at each step
// (cannot make assumptions when cycling over beads!)
neighbor->delay = 0;
neighbor->every = 1;
}
// clang-format off
void FixIPI::initial_integrate(int /*vflag*/)
{
/* This is called at the beginning of the integration loop,

3
src/REPLICA/fix_pimd.cpp
View File

@ -102,6 +102,9 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd", arg[i]));
}
if (strcmp(update->unit_style, "lj") == 0)
error->all(FLERR, "Fix pimd does not support lj units");
/* Initiation */
size_peratom_cols = 12 * nhc_nchain + 3;

184
src/REPLICA/neb.cpp
View File

@ -108,9 +108,9 @@ NEB::~NEB()
void NEB::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"NEB command before simulation box is defined");
error->universe_all(FLERR,"NEB command before simulation box is defined");
if (narg < 6) error->universe_all(FLERR,"Illegal NEB command");
if (narg < 6) error->universe_all(FLERR,"Illegal NEB command: missing argument(s)");
etol = utils::numeric(FLERR,arg[0],false,lmp);
ftol = utils::numeric(FLERR,arg[1],false,lmp);
@ -120,11 +120,14 @@ void NEB::command(int narg, char **arg)
// error checks
if (etol < 0.0) error->all(FLERR,"Illegal NEB command");
if (ftol < 0.0) error->all(FLERR,"Illegal NEB command");
if (nevery <= 0) error->universe_all(FLERR,"Illegal NEB command");
if (n1steps % nevery || n2steps % nevery)
error->universe_all(FLERR,"Illegal NEB command");
if (etol < 0.0) error->universe_all(FLERR, fmt::format("Illegal NEB energy tolerance: {}", etol));
if (ftol < 0.0) error->universe_all(FLERR, fmt::format("Illegal NEB force tolerance: {}", ftol));
if (nevery <= 0)
error->universe_all(FLERR,fmt::format("Illegal NEB command every parameter: {}", nevery));
if (n1steps % nevery)
error->all(FLERR, fmt::format("NEB N1 value {} incompatible with every {}", n1steps, nevery));
if (n2steps % nevery)
error->all(FLERR, fmt::format("NEB N2 value {} incompatible with every {}", n2steps, nevery));
// replica info
@ -136,26 +139,38 @@ void NEB::command(int narg, char **arg)
// error checks
if (nreplica == 1) error->all(FLERR,"Cannot use NEB with a single replica");
if (nreplica == 1) error->universe_all(FLERR,"Cannot use NEB with a single replica");
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"Cannot use NEB unless atom map exists");
error->universe_all(FLERR,"Cannot use NEB without an atom map");
// process file-style setting to setup initial configs for all replicas
if (strcmp(arg[5],"final") == 0) {
if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB command");
inpfile = arg[6];
readfile(inpfile,0);
} else if (strcmp(arg[5],"each") == 0) {
if (narg != 7 && narg !=8) error->universe_all(FLERR,"Illegal NEB command");
inpfile = arg[6];
readfile(inpfile,1);
} else if (strcmp(arg[5],"none") == 0) {
if (narg != 6 && narg !=7) error->universe_all(FLERR,"Illegal NEB command");
} else error->universe_all(FLERR,"Illegal NEB command");
int iarg = 5;
int filecmd = 0;
verbose=false;
if (strcmp(arg[narg-1],"verbose") == 0) verbose=true;
while (iarg < narg) {
if (strcmp(arg[iarg],"final") == 0) {
if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments");
inpfile = arg[iarg+1];
readfile(inpfile,0);
filecmd = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"each") == 0) {
if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments");
inpfile = arg[iarg+1];
readfile(inpfile,1);
filecmd = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"none") == 0) {
filecmd = 1;
++iarg;
} else if (strcmp(arg[iarg],"verbose") == 0) {
verbose=true;
++iarg;
} else error->universe_all(FLERR,fmt::format("Unknown NEB command keyword: {}", arg[iarg]));
}
if (!filecmd) error->universe_all(FLERR, "NEB is missing 'final', 'each', or 'none' keyword");
// run the NEB calculation
run();
@ -176,7 +191,7 @@ void NEB::run()
auto fixes = modify->get_fix_by_style("^neb$");
if (fixes.size() != 1)
error->all(FLERR,"NEB requires use of exactly one fix neb instance");
error->universe_all(FLERR,"NEB requires use of exactly one fix neb instance");
fneb = dynamic_cast<FixNEB *>(fixes[0]);
if (verbose) numall =7;
@ -194,7 +209,7 @@ void NEB::run()
lmp->init();
if (update->minimize->searchflag)
error->all(FLERR,"NEB requires damped dynamics minimizer");
error->universe_all(FLERR,"NEB requires a damped dynamics minimizer");
// setup regular NEB minimization
FILE *uscreen = universe->uscreen;
@ -207,38 +222,43 @@ void NEB::run()
update->endstep = update->laststep = update->firststep + n1steps;
update->nsteps = n1steps;
update->max_eval = n1steps;
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps for NEB");
if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps for NEB");
update->minimize->setup();
if (me_universe == 0) {
if (uscreen) {
fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(uscreen,"\n");
}
if (ulogfile) {
fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(ulogfile,"\n");
}
}
print_status();
@ -292,8 +312,7 @@ void NEB::run()
update->endstep = update->laststep = update->firststep + n2steps;
update->nsteps = n2steps;
update->max_eval = n2steps;
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps");
if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps");
update->minimize->init();
fneb->rclimber = top;
@ -301,34 +320,37 @@ void NEB::run()
if (me_universe == 0) {
if (uscreen) {
fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN "
"pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(uscreen,"\n");
}
if (ulogfile) {
fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ",
"MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT");
for (int i = 1; i <= nreplica; ++i)
fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i));
if (verbose) {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN "
"pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
for (int i = 1; i <= nreplica; ++i) {
auto idx = std::to_string(i);
fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ",
"pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx,
"RepForce"+idx, "MaxAtomForce"+idx);
}
}
fprintf(ulogfile,"\n");
}
}
print_status();
@ -420,7 +442,7 @@ void NEB::readfile(char *file, int flag)
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
if (nlines < 0)
error->all(FLERR,"Incorrectly formatted NEB file");
error->universe_all(FLERR,"Incorrectly formatted NEB file");
}
auto buffer = new char[CHUNK*MAXLINE];
@ -542,11 +564,11 @@ void NEB::open(char *file)
if (platform::has_compress_extension(file)) {
compressed = 1;
fp = platform::compressed_read(file);
if (!fp) error->one(FLERR,"Cannot open compressed file");
if (!fp) error->one(FLERR,"Cannot open compressed file {}: {}", file, utils::getsyserror());
} else fp = fopen(file,"r");
if (fp == nullptr)
error->one(FLERR,"Cannot open file {}: {}",file,utils::getsyserror());
error->one(FLERR,"Cannot open file {}: {}", file, utils::getsyserror());
}
/* ----------------------------------------------------------------------
@ -626,19 +648,19 @@ void NEB::print_status()
if (me_universe == 0) {
constexpr double todeg=180.0/MY_PI;
std::string mesg = fmt::format("{} {:12.8g} {:12.8g} ",update->ntimestep,fmaxreplica,fmaxatom);
mesg += fmt::format("{:12.8g} {:12.8g} {:12.8g} ",gradvnorm0,gradvnorm1,gradvnormc);
mesg += fmt::format("{:12.8g} {:12.8g} {:12.8g} ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:12.8g} {:12.8g} ",rdist[i],all[i][0]);
std::string mesg = fmt::format("{:10} {:<14.8g} {:<14.8g} ",update->ntimestep,fmaxreplica,fmaxatom);
mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",gradvnorm0,gradvnorm1,gradvnormc);
mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:<14.8g} {:<14.8g} ",rdist[i],all[i][0]);
if (verbose) {
mesg += fmt::format("{:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g}",
mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}",
NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg,
all[0][3],freplica[0],fmaxatomInRepl[0]);
for (int i = 1; i < nreplica-1; i++)
mesg += fmt::format("{:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g}",
mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}",
180-acos(all[i][4])*todeg,180-acos(all[i][5])*todeg,
180-acos(all[i][6])*todeg,all[i][3],freplica[i],fmaxatomInRepl[i]);
mesg += fmt::format("{:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g} {:12.5g}",
mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}",
NAN,180-acos(all[nreplica-1][5])*todeg,NAN,all[nreplica-1][3],
freplica[nreplica-1],fmaxatomInRepl[nreplica-1]);
}
@ -650,4 +672,8 @@ void NEB::print_status()
fflush(universe->ulogfile);
}
}
if (verbose) {
delete[] freplica;
delete[] fmaxatomInRepl;
}
}

10
src/angle_hybrid.cpp
View File

@ -306,6 +306,16 @@ void AngleHybrid::coeff(int narg, char **arg)
void AngleHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->nangletypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Angle hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

10
src/bond_hybrid.cpp
View File

@ -334,6 +334,16 @@ void BondHybrid::coeff(int narg, char **arg)
void BondHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->nbondtypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Bond hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();

10
src/dihedral_hybrid.cpp
View File

@ -313,6 +313,16 @@ void DihedralHybrid::coeff(int narg, char **arg)
void DihedralHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->ndihedraltypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Dihedral hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

10
src/improper_hybrid.cpp
View File

@ -305,6 +305,16 @@ void ImproperHybrid::coeff(int narg, char **arg)
void ImproperHybrid::init_style()
{
// error if sub-style is not used
int used;
for (int istyle = 0; istyle < nstyles; ++istyle) {
used = 0;
for (int itype = 1; itype <= atom->nimpropertypes; ++itype)
if (map[itype] == istyle) used = 1;
if (used == 0) error->all(FLERR, "Improper hybrid sub-style {} is not used", keywords[istyle]);
}
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

2
src/input.cpp
View File

@ -1724,7 +1724,7 @@ void Input::pair_coeff()
if (force->pair == nullptr) error->all(FLERR,"Pair_coeff command without a pair style");
if (narg < 2) utils::missing_cmd_args(FLERR,"pair_coeff", error);
if (force->pair->one_coeff && ((strcmp(arg[0],"*") != 0) || (strcmp(arg[1],"*") != 0)))
error->all(FLERR,"Pair_coeff must start with * * for this pair style");
error->all(FLERR,"Pair_coeff must start with * * for pair style {}", force->pair_style);
char *newarg0 = utils::expand_type(FLERR, arg[0], Atom::ATOM, lmp);
if (newarg0) arg[0] = newarg0;

2
src/pair.cpp
View File

@ -821,7 +821,7 @@ void Pair::map_element2type(int narg, char **arg, bool update_setflag)
// elements = list of element names
if (narg != ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR,"Number of element to type mappings does not match number of atom types");
if (elements) {
for (i = 0; i < nelements; i++) delete[] elements[i];

19
src/pair_hybrid.cpp
View File

@ -271,10 +271,9 @@ void PairHybrid::allocate()
void PairHybrid::settings(int narg, char **arg)
{
if (narg < 1) error->all(FLERR,"Illegal pair_style command");
if (narg < 1) utils::missing_cmd_args(FLERR, "pair_style hybrid", error);
if (lmp->kokkos && !utils::strmatch(force->pair_style,"^hybrid.*/kk$"))
error->all(FLERR,"Must use pair_style {}/kk with Kokkos",
force->pair_style);
error->all(FLERR,"Must use pair_style {}/kk with Kokkos", force->pair_style);
// delete old lists, since cannot just change settings
@ -326,9 +325,9 @@ void PairHybrid::settings(int narg, char **arg)
nstyles = 0;
while (iarg < narg) {
if (utils::strmatch(arg[iarg],"^hybrid"))
error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument");
error->all(FLERR,"Pair style hybrid cannot have hybrid as a sub-style");
if (strcmp(arg[iarg],"none") == 0)
error->all(FLERR,"Pair style hybrid cannot have none as an argument");
error->all(FLERR,"Pair style hybrid cannot have none as a sub-style");
styles[nstyles] = force->new_pair(arg[iarg],1,dummy);
keywords[nstyles] = force->store_style(arg[iarg],0);
@ -478,7 +477,7 @@ void PairHybrid::init_svector()
void PairHybrid::coeff(int narg, char **arg)
{
if (narg < 3) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 3) utils::missing_cmd_args(FLERR,"pair_coeff", error);
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -497,7 +496,7 @@ void PairHybrid::coeff(int narg, char **arg)
if (strcmp(arg[2],keywords[m]) == 0) {
if (multiple[m]) {
multflag = 1;
if (narg < 4) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 4) utils::missing_cmd_args(FLERR, "pair_coeff", error);
if (multiple[m] == utils::inumeric(FLERR,arg[3],false,lmp)) break;
else continue;
} else break;
@ -507,7 +506,7 @@ void PairHybrid::coeff(int narg, char **arg)
int none = 0;
if (m == nstyles) {
if (strcmp(arg[2],"none") == 0) none = 1;
else error->all(FLERR,"Pair coeff for hybrid has invalid style: {}",arg[2]);
else error->all(FLERR,"Pair coeff for hybrid has invalid style: {}", arg[2]);
}
// move 1st/2nd args to 2nd/3rd args
@ -527,7 +526,7 @@ void PairHybrid::coeff(int narg, char **arg)
if (!none && styles[m]->one_coeff) {
if ((strcmp(arg[0],"*") != 0) || (strcmp(arg[1],"*") != 0))
error->all(FLERR,"Incorrect args for pair coefficients");
error->all(FLERR,"Pair_coeff must start with * * for sub-style {}", keywords[m]);
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
if (nmap[i][j] && map[i][j][0] == m) {
@ -578,7 +577,7 @@ void PairHybrid::init_style()
for (jtype = itype; jtype <= ntypes; jtype++)
for (m = 0; m < nmap[itype][jtype]; m++)
if (map[itype][jtype][m] == istyle) used = 1;
if (used == 0) error->all(FLERR,"Pair hybrid sub-style is not used");
if (used == 0) error->all(FLERR,"Pair hybrid sub-style {} is not used", keywords[istyle]);
}
// The GPU library uses global data for each pair style, so the