Compare commits

..

39 Commits

Author SHA1 Message Date
6fd8b2b177 Merge pull request #3122 from akohlmey/maintenance-2021-09-29
Third round of maintenance fixes for the stable release
2022-03-24 14:20:52 -04:00
6edaf42b3d fix temperature initialization bug in KOKKOS nose-hoover code 2022-03-24 11:44:24 -04:00
79c047487d fix parallel execution bug for shell command 2022-03-24 07:38:44 -04:00
ac5acb9abf update threebody example 2022-03-24 07:31:02 -04:00
87fbbd3b13 small kokkos fixes from upstream 2022-03-24 07:18:24 -04:00
8ac0ec6473 Changes needed to compile LAMMPS with latest Kokkos develop 2022-03-24 06:09:03 -04:00
8acba74c4d correct input to load potential file from local folder 2022-03-22 22:32:39 -04:00
34bcbdf41d update extep potential file 2022-03-22 22:31:48 -04:00
d519ca0213 add missing reaxff files to purge list 2022-03-21 14:34:14 -04:00
a392e8dc09 accept infile with 0 lines, so we can create a template from the restart 2022-03-21 00:33:40 -04:00
a4d4f77bc2 run setup_bodies_dynamic() before processing infile in case that is not resetting all data 2022-03-21 00:32:49 -04:00
83a8f72d83 fix off-by-one bug when writing restart files for rigid bodies 2022-03-20 19:14:13 -04:00
3c54b56cfe update overlooked date stamp 2022-03-19 21:00:14 -04:00
ff1a08f148 fixes to CMake build for ML-QUIP package from upstream 2022-03-17 18:07:12 -04:00
5a53b0fc03 import python3 compatibility changes to tools/python from upstream 2022-03-16 13:24:53 -04:00
e550600ebe Error fixed. Epsilon and sigma must also be symmetric 2022-03-16 09:09:52 -04:00
7cb13be52a fix bug where it was not possible to use an absolute path for write_coeff 2022-03-16 09:08:47 -04:00
ab56d7ecd7 augment cmake library search path to include the CUDA stubs library folder
this will help configuring and compiling LAMMPS with CUDA support on
machines where there is no CUDA driver installed
2022-03-10 23:02:57 -05:00
bd6ac3ee6d for 2d systems, rigid bodies always have a moment of inertia and no DOFs need to be subtracted 2022-03-02 16:41:35 -05:00
27ca0a8f41 trigger building an "intel" style neighbor list so that buffers are allocated 2022-02-27 14:50:48 -05:00
f688b9b6b5 use consistent names, avoid memory leaks, fix off-by-1 error in fourier dihedral 2022-02-27 12:25:32 -05:00
16c61b3cc0 add support for plumed 2.6.5, 2.6.6, 2.7.3, 2.7.4, and 2.8.0 (default 2.7.4) 2022-02-25 16:37:00 -05:00
fb480f22fc make cythonize detection compatible with /bin/dash on ubunutu 2022-02-24 21:24:04 -05:00
d0507559a4 when updating ML-IAP due to adding/removing PYTHON we need to delete and re-add cythonize support 2022-02-24 20:40:55 -05:00
ali
58eb331b08 Python 3 compatibility for log commands in tools/python 2022-02-23 10:22:29 -05:00
c68015ca87 Bug fix for Intel package skip lists with multiple runs. 2022-02-18 05:11:34 -05:00
583c22d6e0 update tools/eam_database from upstream 2022-02-16 11:46:11 -05:00
58a4694d92 Remove incorrect error check in ReaxFF 2022-02-11 16:19:00 -05:00
97cf345528 don't allow exceptions to "escape" a destructor 2022-02-10 21:13:26 -05:00
0658abbdd4 silence possible warnings about missing files on "make clean-all" 2022-02-10 21:10:34 -05:00
72026a58bf make certain that "offset" is always initialized 2022-02-10 21:05:12 -05:00
7152231a10 plug memory leak 2022-02-10 20:56:51 -05:00
8fe8a667b6 update create.f with changes from NIST database
also add parameters for Cr and document in README file and change
the code to create output files with .eam.alloy extension
2022-02-10 20:45:16 -05:00
560c543e69 add extra communication of special neighbors when using angle constraints 2022-02-10 20:44:39 -05:00
c5e6650924 import bugfixes for crashes and memory leaks in MSM kspace style from develop 2022-02-10 20:36:35 -05:00
10373ea5c9 avoid failures with "most" presets 2022-02-10 20:11:00 -05:00
992b1cf582 label as update #3 2022-01-25 07:42:00 -05:00
1505f3de06 fix tag caching issue in INTEL package 2022-01-25 07:41:37 -05:00
566efe04f2 always fall back to using the .so extension if available in the LAMMPS module folder 2022-01-19 10:12:50 -05:00
84 changed files with 2075 additions and 918 deletions

View File

@ -1,10 +1,11 @@
find_package(ZLIB REQUIRED)
target_link_libraries(lammps PRIVATE ZLIB::ZLIB)
find_package(PkgConfig REQUIRED)
pkg_check_modules(Zstd IMPORTED_TARGET libzstd>=1.4)
if(Zstd_FOUND)
find_package(PkgConfig QUIET)
if(PkgConfig_FOUND)
pkg_check_modules(Zstd IMPORTED_TARGET libzstd>=1.4)
if(Zstd_FOUND)
target_compile_definitions(lammps PRIVATE -DLAMMPS_ZSTD)
target_link_libraries(lammps PRIVATE PkgConfig::Zstd)
endif()
endif()

View File

@ -30,7 +30,15 @@ file(GLOB GPU_LIB_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/[^.]*.cpp)
file(MAKE_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/gpu)
if(GPU_API STREQUAL "CUDA")
find_package(CUDA REQUIRED)
find_package(CUDA QUIET)
# augment search path for CUDA toolkit libraries to include the stub versions. Needed to find libcuda.so on machines without a CUDA driver installation
if(CUDA_FOUND)
set(CMAKE_LIBRARY_PATH "${CUDA_TOOLKIT_ROOT_DIR}/lib64/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib64;${CUDA_TOOLKIT_ROOT_DIR}/lib;${CMAKE_LIBRARY_PATH}")
find_package(CUDA REQUIRED)
else()
message(FATAL_ERROR "CUDA Toolkit not found")
endif()
find_program(BIN2C bin2c)
if(NOT BIN2C)
message(FATAL_ERROR "Could not find bin2c, use -DBIN2C=/path/to/bin2c to help cmake finding it.")

View File

@ -32,7 +32,8 @@ if(DOWNLOAD_QUIP)
foreach(flag ${LAPACK_LIBRARIES})
set(temp "${temp} ${flag}")
endforeach()
set(temp "${temp}\n")
# Fix cmake crashing when MATH_LINKOPTS not set, required for e.g. recent Cray Programming Environment
set(temp "${temp} -L/_DUMMY_PATH_\n")
set(temp "${temp}PYTHON=python\nPIP=pip\nEXTRA_LINKOPTS=\n")
set(temp "${temp}HAVE_CP2K=0\nHAVE_VASP=0\nHAVE_TB=0\nHAVE_PRECON=1\nHAVE_LOTF=0\nHAVE_ONIOM=0\n")
set(temp "${temp}HAVE_LOCAL_E_MIX=0\nHAVE_QC=0\nHAVE_GAP=1\nHAVE_DESCRIPTORS_NONCOMMERCIAL=1\n")
@ -50,6 +51,7 @@ if(DOWNLOAD_QUIP)
GIT_TAG origin/public
GIT_SHALLOW YES
GIT_PROGRESS YES
GIT_SUBMODULES "src/fox;src/GAP"
PATCH_COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_BINARY_DIR}/quip.config <SOURCE_DIR>/arch/Makefile.lammps
CONFIGURE_COMMAND env QUIP_ARCH=lammps make config
BUILD_COMMAND env QUIP_ARCH=lammps make libquip

View File

@ -54,8 +54,8 @@ if(DOWNLOAD_PLUMED)
set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/libplumedWrapper.a")
endif()
set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.7.2/plumed-src-2.7.2.tgz" CACHE STRING "URL for PLUMED tarball")
set(PLUMED_MD5 "cfa0b4dd90a81c25d3302e8d97bfeaea" CACHE STRING "MD5 checksum of PLUMED tarball")
set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.7.4/plumed-src-2.7.4.tgz" CACHE STRING "URL for PLUMED tarball")
set(PLUMED_MD5 "858e0b6aed173748fc85b6bc8a9dcb3e" CACHE STRING "MD5 checksum of PLUMED tarball")
mark_as_advanced(PLUMED_URL)
mark_as_advanced(PLUMED_MD5)

View File

@ -48,7 +48,6 @@ set(ALL_PACKAGES
PHONON
PLUGIN
POEMS
PYTHON
QEQ
REACTION
REAXFF

View File

@ -278,16 +278,20 @@ eam database tool
-----------------------------
The tools/eam_database directory contains a Fortran program that will
generate EAM alloy setfl potential files for any combination of 16
generate EAM alloy setfl potential files for any combination of 17
elements: Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti,
Zr. The files can then be used with the :doc:`pair_style eam/alloy <pair_eam>` command.
Zr, Cr. The files can then be used with the :doc:`pair_style eam/alloy <pair_eam>` command.
The tool is authored by Xiaowang Zhou (Sandia), xzhou at sandia.gov,
and is based on his paper:
with updates from Lucas Hale (NIST) lucas.hale at nist.gov and is based on his paper:
X. W. Zhou, R. A. Johnson, and H. N. G. Wadley, Phys. Rev. B, 69,
144113 (2004).
The parameters for Cr were taken from:
Lin Z B, Johnson R A and Zhigilei L V, Phys. Rev. B 77 214108 (2008).
----------
.. _eamgn:

View File

@ -0,0 +1 @@
../../../potentials/BN.extep

View File

@ -15,7 +15,7 @@ neigh_modify check yes
# Potential
pair_style extep
pair_coeff * * ../../../../potentials/BN.extep B N
pair_coeff * * BN.extep B N
# Output
thermo 10

View File

@ -0,0 +1 @@
../../potentials/SiC.tersoff.zbl

View File

@ -7,7 +7,7 @@ units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
atom_modify sort 0 0.0
# temperature
@ -35,23 +35,23 @@ region myreg block 0 4 &
create_box 8 myreg
create_atoms 1 region myreg &
basis 1 1 &
basis 2 2 &
basis 3 3 &
basis 4 4 &
basis 5 5 &
basis 6 6 &
basis 7 7 &
basis 8 8
basis 1 1 &
basis 2 2 &
basis 3 3 &
basis 4 4 &
basis 5 5 &
basis 6 6 &
basis 7 7 &
basis 8 8
mass * 28.06
velocity all create $t 5287287 loop geom
velocity all create $t 5287287 loop geom
# Equilibrate using Stillinger-Weber model for silicon
pair_style sw
pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
thermo_style custom step temp epair etotal econserve press
thermo 10
@ -61,15 +61,15 @@ neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
run 100
write_restart restart.equil
write_restart restart.equil
# Test Stillinger-Weber model for Cd/Te/Zn/Se/Hg/S
clear
read_restart restart.equil
read_restart restart.equil
pair_style sw
pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
thermo_style custom step temp epair etotal econserve press
thermo 10
@ -82,10 +82,10 @@ run 100
# Test Vashishta model for In/P
clear
read_restart restart.equil
read_restart restart.equil
pair_style vashishta
pair_coeff * * InP.vashishta In In In In P P P P
pair_coeff * * InP.vashishta In In In In P P P P
thermo_style custom step temp epair etotal econserve press
thermo 10
@ -98,13 +98,13 @@ run 100
# Test Tersoff model for B/N/C
clear
read_restart restart.equil
read_restart restart.equil
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
pair_style tersoff
pair_coeff * * BNC.tersoff N N N C B B C B
pair_coeff * * BNC.tersoff N N N C B B C B
thermo_style custom step temp epair etotal econserve press
thermo 10
@ -114,3 +114,23 @@ neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
run 100
# Test Tersoff/ZBL model for SiC
clear
read_restart restart.equil
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
pair_style tersoff/zbl
pair_coeff * * SiC.tersoff.zbl C C C C Si Si Si Si
thermo_style custom step temp epair etotal econserve press
thermo 10
fix 1 all nvt temp $t $t 0.1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
run 100
shell rm restart.equil

View File

@ -1,4 +1,4 @@
LAMMPS (24 Dec 2020)
LAMMPS (29 Sep 2021 - Update 3)
using 1 OpenMP thread(s) per MPI task
# Simple regression tests for threebody potentials
@ -9,7 +9,7 @@ units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
atom_modify sort 0 0.0
# temperature
@ -27,19 +27,20 @@ region myreg block 0 4 0 4
create_box 8 myreg
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 1 by 1 MPI processor grid
create_atoms 1 region myreg basis 1 1 basis 2 2 basis 3 3 basis 4 4 basis 5 5 basis 6 6 basis 7 7 basis 8 8
create_atoms 1 region myreg basis 1 1 basis 2 2 basis 3 3 basis 4 4 basis 5 5 basis 6 6 basis 7 7 basis 8 8
Created 512 atoms
create_atoms CPU = 0.001 seconds
using lattice units in orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
create_atoms CPU = 0.000 seconds
mass * 28.06
velocity all create $t 5287287 loop geom
velocity all create 1800 5287287 loop geom
velocity all create $t 5287287 loop geom
velocity all create 1800 5287287 loop geom
# Equilibrate using Stillinger-Weber model for silicon
pair_style sw
pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
Reading sw potential file Si.sw with DATE: 2007-06-11
thermo_style custom step temp epair etotal econserve press
@ -76,20 +77,20 @@ Step Temp E_pair TotEng Econserve Press
80 800.80221 -2146.1371 -2093.2426 -2101.313 11995.66
90 1293.9689 -2176.9021 -2091.4329 -2101.3848 11692.45
100 1112.9699 -2162.7259 -2089.2121 -2101.3478 12263.758
Loop time of 0.157871 on 1 procs for 100 steps with 512 atoms
Loop time of 0.093281 on 1 procs for 100 steps with 512 atoms
Performance: 54.728 ns/day, 0.439 hours/ns, 633.430 timesteps/s
Performance: 92.623 ns/day, 0.259 hours/ns, 1072.029 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 | 0.14704 | 0.14704 | 0.14704 | 0.0 | 93.14
Neigh | 0.00247 | 0.00247 | 0.00247 | 0.0 | 1.56
Comm | 0.0024729 | 0.0024729 | 0.0024729 | 0.0 | 1.57
Output | 0.0002656 | 0.0002656 | 0.0002656 | 0.0 | 0.17
Modify | 0.0050237 | 0.0050237 | 0.0050237 | 0.0 | 3.18
Other | | 0.0006011 | | | 0.38
Pair | 0.090256 | 0.090256 | 0.090256 | 0.0 | 96.76
Neigh | 0.0015078 | 0.0015078 | 0.0015078 | 0.0 | 1.62
Comm | 0.00045896 | 0.00045896 | 0.00045896 | 0.0 | 0.49
Output | 8.3447e-05 | 8.3447e-05 | 8.3447e-05 | 0.0 | 0.09
Modify | 0.00072384 | 0.00072384 | 0.00072384 | 0.0 | 0.78
Other | | 0.0002506 | | | 0.27
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -105,25 +106,25 @@ Ave neighs/atom = 27.320312
Neighbor list builds = 2
Dangerous builds = 0
write_restart restart.equil
write_restart restart.equil
System init for write_restart ...
# Test Stillinger-Weber model for Cd/Te/Zn/Se/Hg/S
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
read_restart restart.equil
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.001 seconds
read_restart CPU = 0.000 seconds
pair_style sw
pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
Reading sw potential file CdTeZnSeHgS0.sw with DATE: 2013-08-09
thermo_style custom step temp epair etotal econserve press
@ -163,20 +164,20 @@ Step Temp E_pair TotEng Econserve Press
180 1856.1197 -657.14338 -534.54309 -564.48754 488372.27
190 1346.1107 -621.42431 -532.5111 -564.38065 511750.04
200 1919.5266 -657.26587 -530.47743 -564.47797 488684.56
Loop time of 0.455825 on 1 procs for 100 steps with 512 atoms
Loop time of 0.245572 on 1 procs for 100 steps with 512 atoms
Performance: 18.955 ns/day, 1.266 hours/ns, 219.382 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
Performance: 35.183 ns/day, 0.682 hours/ns, 407.212 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 | 0.44091 | 0.44091 | 0.44091 | 0.0 | 96.73
Neigh | 0.0054555 | 0.0054555 | 0.0054555 | 0.0 | 1.20
Comm | 0.0035784 | 0.0035784 | 0.0035784 | 0.0 | 0.79
Output | 0.00024486 | 0.00024486 | 0.00024486 | 0.0 | 0.05
Modify | 0.0050471 | 0.0050471 | 0.0050471 | 0.0 | 1.11
Other | | 0.000592 | | | 0.13
Pair | 0.24139 | 0.24139 | 0.24139 | 0.0 | 98.30
Neigh | 0.0027068 | 0.0027068 | 0.0027068 | 0.0 | 1.10
Comm | 0.00051188 | 0.00051188 | 0.00051188 | 0.0 | 0.21
Output | 0.00010395 | 0.00010395 | 0.00010395 | 0.0 | 0.04
Modify | 0.00059605 | 0.00059605 | 0.00059605 | 0.0 | 0.24
Other | | 0.0002608 | | | 0.11
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -196,18 +197,18 @@ Dangerous builds = 0
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
read_restart restart.equil
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.001 seconds
read_restart CPU = 0.000 seconds
pair_style vashishta
pair_coeff * * InP.vashishta In In In In P P P P
pair_coeff * * InP.vashishta In In In In P P P P
Reading vashishta potential file InP.vashishta with DATE: 2015-10-14
thermo_style custom step temp epair etotal econserve press
@ -247,20 +248,20 @@ Step Temp E_pair TotEng Econserve Press
180 1302.9041 -1491.7765 -1405.7172 -1435.8971 249514.04
190 1332.3326 -1491.5271 -1403.524 -1435.9213 227537.99
200 1352.1813 -1490.4513 -1401.1371 -1435.9049 207626.42
Loop time of 0.217808 on 1 procs for 100 steps with 512 atoms
Loop time of 0.111899 on 1 procs for 100 steps with 512 atoms
Performance: 39.668 ns/day, 0.605 hours/ns, 459.121 timesteps/s
98.2% CPU use with 1 MPI tasks x 1 OpenMP threads
Performance: 77.212 ns/day, 0.311 hours/ns, 893.662 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.19635 | 0.19635 | 0.19635 | 0.0 | 90.15
Neigh | 0.01054 | 0.01054 | 0.01054 | 0.0 | 4.84
Comm | 0.0051923 | 0.0051923 | 0.0051923 | 0.0 | 2.38
Output | 0.00027919 | 0.00027919 | 0.00027919 | 0.0 | 0.13
Modify | 0.0048637 | 0.0048637 | 0.0048637 | 0.0 | 2.23
Other | | 0.0005858 | | | 0.27
Pair | 0.10539 | 0.10539 | 0.10539 | 0.0 | 94.18
Neigh | 0.0049229 | 0.0049229 | 0.0049229 | 0.0 | 4.40
Comm | 0.00068307 | 0.00068307 | 0.00068307 | 0.0 | 0.61
Output | 6.1989e-05 | 6.1989e-05 | 6.1989e-05 | 0.0 | 0.06
Modify | 0.00058532 | 0.00058532 | 0.00058532 | 0.0 | 0.52
Other | | 0.0002604 | | | 0.23
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -280,28 +281,28 @@ Dangerous builds = 0
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
read_restart restart.equil
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.001 seconds
read_restart CPU = 0.000 seconds
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap
Changing box ...
orthogonal box = (4.3448000 0.0000000 0.0000000) to (17.379200 21.724000 21.724000)
orthogonal box = (4.3448000 4.3448000 0.0000000) to (17.379200 17.379200 21.724000)
orthogonal box = (4.3448000 4.3448000 4.3448000) to (17.379200 17.379200 17.379200)
pair_style tersoff
pair_coeff * * BNC.tersoff N N N C B B C B
pair_coeff * * BNC.tersoff N N N C B B C B
Reading tersoff potential file BNC.tersoff with DATE: 2013-03-21
thermo_style custom step temp epair etotal econserve press
@ -341,20 +342,20 @@ Step Temp E_pair TotEng Econserve Press
180 1337.4358 -3254.9844 -3166.6442 -3196.8222 1880420.9
190 1441.8052 -3259.0364 -3163.8023 -3196.3556 1904512.1
200 1569.0317 -3265.0089 -3161.3714 -3196.3328 1899462.7
Loop time of 0.487425 on 1 procs for 100 steps with 512 atoms
Loop time of 0.097734 on 1 procs for 100 steps with 512 atoms
Performance: 17.726 ns/day, 1.354 hours/ns, 205.160 timesteps/s
99.1% CPU use with 1 MPI tasks x 1 OpenMP threads
Performance: 88.403 ns/day, 0.271 hours/ns, 1023.186 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.47762 | 0.47762 | 0.47762 | 0.0 | 97.99
Neigh | 0.0014286 | 0.0014286 | 0.0014286 | 0.0 | 0.29
Comm | 0.0024068 | 0.0024068 | 0.0024068 | 0.0 | 0.49
Output | 0.00028992 | 0.00028992 | 0.00028992 | 0.0 | 0.06
Modify | 0.0050635 | 0.0050635 | 0.0050635 | 0.0 | 1.04
Other | | 0.0006182 | | | 0.13
Pair | 0.095481 | 0.095481 | 0.095481 | 0.0 | 97.69
Neigh | 0.000772 | 0.000772 | 0.000772 | 0.0 | 0.79
Comm | 0.00046158 | 0.00046158 | 0.00046158 | 0.0 | 0.47
Output | 6.7949e-05 | 6.7949e-05 | 6.7949e-05 | 0.0 | 0.07
Modify | 0.00068784 | 0.00068784 | 0.00068784 | 0.0 | 0.70
Other | | 0.0002635 | | | 0.27
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -370,4 +371,99 @@ Ave neighs/atom = 28.523438
Neighbor list builds = 1
Dangerous builds = 0
# Test Tersoff/ZBL model for SiC
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
Reading restart file ...
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.000 seconds
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap
Changing box ...
orthogonal box = (4.3448000 0.0000000 0.0000000) to (17.379200 21.724000 21.724000)
orthogonal box = (4.3448000 4.3448000 0.0000000) to (17.379200 17.379200 21.724000)
orthogonal box = (4.3448000 4.3448000 4.3448000) to (17.379200 17.379200 17.379200)
pair_style tersoff/zbl
pair_coeff * * SiC.tersoff.zbl C C C C Si Si Si Si
Reading tersoff/zbl potential file SiC.tersoff.zbl with DATE: 2009-04-15
thermo_style custom step temp epair etotal econserve press
thermo 10
fix 1 all nvt temp $t $t 0.1
fix 1 all nvt temp 1800 $t 0.1
fix 1 all nvt temp 1800 1800 0.1
Resetting global fix info from restart file:
fix style: nvt, fix ID: 1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
run 100
All restart file global fix info was re-assigned
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4
ghost atom cutoff = 4
binsize = 2, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair tersoff/zbl, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.002 | 3.002 | 3.002 Mbytes
Step Temp E_pair TotEng Econserve Press
100 1112.9699 7067.9634 7141.4772 7129.3415 17683957
110 1676.669 7033.1458 7143.893 7128.6921 17837566
120 2450.2667 6982.2491 7144.094 7126.9524 18220027
130 2726.9659 6964.1219 7144.2432 7126.7678 18230324
140 2729.421 6962.7393 7143.0228 7127.2074 18176317
150 2738.5449 6959.1761 7140.0623 7127.6671 18068370
160 2687.2419 6958.1183 7135.6158 7127.8492 18156214
170 2697.7325 6952.1482 7130.3387 7127.7898 17978251
180 2577.9885 6954.5611 7124.8422 7127.5615 18068920
190 2502.6928 6954.4558 7119.7635 7127.67 18049652
200 2517.4866 6947.962 7114.2469 7127.1972 18209451
Loop time of 0.783169 on 1 procs for 100 steps with 512 atoms
Performance: 11.032 ns/day, 2.175 hours/ns, 127.686 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.78056 | 0.78056 | 0.78056 | 0.0 | 99.67
Neigh | 0.0011299 | 0.0011299 | 0.0011299 | 0.0 | 0.14
Comm | 0.00051332 | 0.00051332 | 0.00051332 | 0.0 | 0.07
Output | 9.2268e-05 | 9.2268e-05 | 9.2268e-05 | 0.0 | 0.01
Modify | 0.00060058 | 0.00060058 | 0.00060058 | 0.0 | 0.08
Other | | 0.0002706 | | | 0.03
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1542.00 ave 1542 max 1542 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 30142.0 ave 30142 max 30142 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 30142
Ave neighs/atom = 58.871094
Neighbor list builds = 1
Dangerous builds = 0
shell rm restart.equil
Total wall time: 0:00:01

View File

@ -1,4 +1,4 @@
LAMMPS (24 Dec 2020)
LAMMPS (29 Sep 2021 - Update 3)
using 1 OpenMP thread(s) per MPI task
# Simple regression tests for threebody potentials
@ -9,7 +9,7 @@ units metal
atom_style atomic
atom_modify map array
boundary p p p
atom_modify sort 0 0.0
atom_modify sort 0 0.0
# temperature
@ -26,20 +26,21 @@ region myreg block 0 4 0 4
create_box 8 myreg
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 2 by 2 MPI processor grid
create_atoms 1 region myreg basis 1 1 basis 2 2 basis 3 3 basis 4 4 basis 5 5 basis 6 6 basis 7 7 basis 8 8
1 by 1 by 1 MPI processor grid
create_atoms 1 region myreg basis 1 1 basis 2 2 basis 3 3 basis 4 4 basis 5 5 basis 6 6 basis 7 7 basis 8 8
Created 512 atoms
create_atoms CPU = 0.074 seconds
using lattice units in orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
create_atoms CPU = 0.000 seconds
mass * 28.06
velocity all create $t 5287287 loop geom
velocity all create 1800 5287287 loop geom
velocity all create $t 5287287 loop geom
velocity all create 1800 5287287 loop geom
# Equilibrate using Stillinger-Weber model for silicon
pair_style sw
pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
Reading sw potential file Si.sw with DATE: 2007-06-11
thermo_style custom step temp epair etotal econserve press
@ -63,7 +64,7 @@ Neighbor list info ...
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.958 | 2.958 | 2.958 Mbytes
Per MPI rank memory allocation (min/avg/max) = 2.983 | 2.983 | 2.983 Mbytes
Step Temp E_pair TotEng Econserve Press
0 1800 -2220.3392 -2101.4457 -2101.4457 12358.626
10 1006.0192 -2167.7053 -2101.2558 -2101.3286 13892.426
@ -76,54 +77,54 @@ Step Temp E_pair TotEng Econserve Press
80 800.80221 -2146.1371 -2093.2426 -2101.313 11995.66
90 1293.9689 -2176.9021 -2091.4329 -2101.3848 11692.45
100 1112.9699 -2162.7259 -2089.2121 -2101.3478 12263.758
Loop time of 0.0998364 on 4 procs for 100 steps with 512 atoms
Loop time of 0.089642 on 1 procs for 100 steps with 512 atoms
Performance: 86.542 ns/day, 0.277 hours/ns, 1001.639 timesteps/s
81.4% CPU use with 4 MPI tasks x 1 OpenMP threads
Performance: 96.383 ns/day, 0.249 hours/ns, 1115.548 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.037337 | 0.049389 | 0.069239 | 5.9 | 49.47
Neigh | 0.00067854 | 0.00068814 | 0.00070286 | 0.0 | 0.69
Comm | 0.025239 | 0.04504 | 0.056869 | 6.1 | 45.11
Output | 0.00015712 | 0.00082219 | 0.0028148 | 0.0 | 0.82
Modify | 0.0014369 | 0.0015754 | 0.0016632 | 0.2 | 1.58
Other | | 0.002321 | | | 2.33
Pair | 0.086619 | 0.086619 | 0.086619 | 0.0 | 96.63
Neigh | 0.0015211 | 0.0015211 | 0.0015211 | 0.0 | 1.70
Comm | 0.000458 | 0.000458 | 0.000458 | 0.0 | 0.51
Output | 7.987e-05 | 7.987e-05 | 7.987e-05 | 0.0 | 0.09
Modify | 0.00073361 | 0.00073361 | 0.00073361 | 0.0 | 0.82
Other | | 0.0002301 | | | 0.26
Nlocal: 128.000 ave 132 max 125 min
Histogram: 1 1 0 0 0 1 0 0 0 1
Nghost: 525.000 ave 528 max 521 min
Histogram: 1 0 0 0 1 0 0 0 1 1
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1017.00 ave 1017 max 1017 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 3497.00 ave 3619 max 3397 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 13988.0 ave 13988 max 13988 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 13988
Ave neighs/atom = 27.320312
Neighbor list builds = 2
Dangerous builds = 0
write_restart restart.equil
write_restart restart.equil
System init for write_restart ...
# Test Stillinger-Weber model for Cd/Te/Zn/Se/Hg/S
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
read_restart restart.equil
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 2 by 2 MPI processor grid
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.001 seconds
read_restart CPU = 0.000 seconds
pair_style sw
pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
Reading sw potential file CdTeZnSeHgS0.sw with DATE: 2013-08-09
thermo_style custom step temp epair etotal econserve press
@ -150,42 +151,42 @@ Neighbor list info ...
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.967 | 2.967 | 2.968 Mbytes
Per MPI rank memory allocation (min/avg/max) = 3.001 | 3.001 | 3.001 Mbytes
Step Temp E_pair TotEng Econserve Press
100 1112.9699 -625.76163 -552.24782 -564.38354 462129.66
100 1112.9699 -625.76163 -552.24781 -564.38354 462129.66
110 1502.8461 -649.55768 -550.29179 -564.45814 463413.45
120 1926.4523 -674.71265 -547.46675 -564.53613 486338.88
130 1152.6663 -621.47265 -545.33681 -564.37203 514892.19
120 1926.4523 -674.71265 -547.46675 -564.53612 486338.88
130 1152.6663 -621.47264 -545.33681 -564.37203 514892.2
140 1762.244 -659.86941 -543.46979 -564.4985 488159.88
150 1767.8665 -657.67179 -540.90079 -564.48386 466721.31
160 1075.2874 -610.1281 -539.10328 -564.36709 470151.9
170 1697.9313 -649.3684 -537.21676 -564.47208 467953.7
180 1856.1197 -657.14338 -534.54309 -564.48754 488372.26
190 1346.1107 -621.42432 -532.5111 -564.38065 511750.03
150 1767.8665 -657.67178 -540.90078 -564.48386 466721.31
160 1075.2874 -610.12809 -539.10328 -564.36709 470151.9
170 1697.9313 -649.3684 -537.21675 -564.47207 467953.71
180 1856.1197 -657.14338 -534.54309 -564.48754 488372.27
190 1346.1107 -621.42431 -532.5111 -564.38065 511750.04
200 1919.5266 -657.26587 -530.47743 -564.47797 488684.56
Loop time of 0.286556 on 4 procs for 100 steps with 512 atoms
Loop time of 0.268183 on 1 procs for 100 steps with 512 atoms
Performance: 30.151 ns/day, 0.796 hours/ns, 348.971 timesteps/s
81.7% CPU use with 4 MPI tasks x 1 OpenMP threads
Performance: 32.217 ns/day, 0.745 hours/ns, 372.880 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 | 0.11093 | 0.139 | 0.16864 | 5.8 | 48.51
Neigh | 0.0014305 | 0.0014756 | 0.0015156 | 0.1 | 0.51
Comm | 0.10154 | 0.12374 | 0.16907 | 7.8 | 43.18
Output | 0.0001862 | 0.00030428 | 0.0006578 | 0.0 | 0.11
Modify | 0.0038164 | 0.019159 | 0.034146 | 10.8 | 6.69
Other | | 0.002872 | | | 1.00
Pair | 0.26374 | 0.26374 | 0.26374 | 0.0 | 98.34
Neigh | 0.0027301 | 0.0027301 | 0.0027301 | 0.0 | 1.02
Comm | 0.00063014 | 0.00063014 | 0.00063014 | 0.0 | 0.23
Output | 8.4639e-05 | 8.4639e-05 | 8.4639e-05 | 0.0 | 0.03
Modify | 0.00072742 | 0.00072742 | 0.00072742 | 0.0 | 0.27
Other | | 0.0002725 | | | 0.10
Nlocal: 128.000 ave 135 max 122 min
Histogram: 1 0 1 0 0 0 1 0 0 1
Nghost: 759.750 ave 770 max 751 min
Histogram: 1 0 0 1 1 0 0 0 0 1
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1428.00 ave 1428 max 1428 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 4336.00 ave 4563 max 4128 min
Histogram: 1 0 1 0 0 0 1 0 0 1
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 17344.0 ave 17344 max 17344 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 17344
Ave neighs/atom = 33.875000
@ -196,18 +197,18 @@ Dangerous builds = 0
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
read_restart restart.equil
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 2 by 2 MPI processor grid
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.001 seconds
read_restart CPU = 0.000 seconds
pair_style vashishta
pair_coeff * * InP.vashishta In In In In P P P P
pair_coeff * * InP.vashishta In In In In P P P P
Reading vashishta potential file InP.vashishta with DATE: 2015-10-14
thermo_style custom step temp epair etotal econserve press
@ -234,7 +235,7 @@ Neighbor list info ...
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.988 | 2.988 | 2.988 Mbytes
Per MPI rank memory allocation (min/avg/max) = 3.025 | 3.025 | 3.025 Mbytes
Step Temp E_pair TotEng Econserve Press
100 1112.9699 -1497.2988 -1423.785 -1435.9207 355619.19
110 1250.545 -1504.5795 -1421.9785 -1435.9786 345188.52
@ -247,29 +248,29 @@ Step Temp E_pair TotEng Econserve Press
180 1302.9041 -1491.7765 -1405.7172 -1435.8971 249514.04
190 1332.3326 -1491.5271 -1403.524 -1435.9213 227537.99
200 1352.1813 -1490.4513 -1401.1371 -1435.9049 207626.42
Loop time of 0.14468 on 4 procs for 100 steps with 512 atoms
Loop time of 0.117875 on 1 procs for 100 steps with 512 atoms
Performance: 59.718 ns/day, 0.402 hours/ns, 691.179 timesteps/s
81.2% CPU use with 4 MPI tasks x 1 OpenMP threads
Performance: 73.298 ns/day, 0.327 hours/ns, 848.357 timesteps/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.047903 | 0.058669 | 0.086091 | 6.6 | 40.55
Neigh | 0.0027876 | 0.002852 | 0.0028808 | 0.1 | 1.97
Comm | 0.034642 | 0.066142 | 0.078599 | 7.1 | 45.72
Output | 0.00018477 | 0.0049147 | 0.019101 | 11.7 | 3.40
Modify | 0.0015709 | 0.0022651 | 0.0029545 | 1.4 | 1.57
Other | | 0.009837 | | | 6.80
Pair | 0.11085 | 0.11085 | 0.11085 | 0.0 | 94.04
Neigh | 0.005235 | 0.005235 | 0.005235 | 0.0 | 4.44
Comm | 0.00077152 | 0.00077152 | 0.00077152 | 0.0 | 0.65
Output | 7.9155e-05 | 7.9155e-05 | 7.9155e-05 | 0.0 | 0.07
Modify | 0.00065637 | 0.00065637 | 0.00065637 | 0.0 | 0.56
Other | | 0.0002811 | | | 0.24
Nlocal: 128.000 ave 131 max 124 min
Histogram: 1 0 0 0 0 1 0 1 0 1
Nghost: 1013.25 ave 1025 max 1002 min
Histogram: 1 1 0 0 0 0 0 0 1 1
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1838.00 ave 1838 max 1838 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 9120.50 ave 9356 max 8868 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 36482.0 ave 36482 max 36482 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 36482
Ave neighs/atom = 71.253906
@ -280,28 +281,28 @@ Dangerous builds = 0
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
read_restart restart.equil
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 2 by 2 MPI processor grid
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.005 seconds
read_restart CPU = 0.000 seconds
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap
Changing box ...
orthogonal box = (4.3448000 0.0000000 0.0000000) to (17.379200 21.724000 21.724000)
orthogonal box = (4.3448000 4.3448000 0.0000000) to (17.379200 17.379200 21.724000)
orthogonal box = (4.3448000 4.3448000 4.3448000) to (17.379200 17.379200 17.379200)
pair_style tersoff
pair_coeff * * BNC.tersoff N N N C B B C B
pair_coeff * * BNC.tersoff N N N C B B C B
Reading tersoff potential file BNC.tersoff with DATE: 2013-03-21
thermo_style custom step temp epair etotal econserve press
@ -328,7 +329,7 @@ Neighbor list info ...
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.948 | 2.948 | 2.948 Mbytes
Per MPI rank memory allocation (min/avg/max) = 2.982 | 2.982 | 2.982 Mbytes
Step Temp E_pair TotEng Econserve Press
100 1112.9699 -3259.7676 -3186.2538 -3198.3895 1912461.3
110 1772.8268 -3301.5479 -3184.4493 -3198.8218 1885295.6
@ -341,33 +342,128 @@ Step Temp E_pair TotEng Econserve Press
180 1337.4358 -3254.9844 -3166.6442 -3196.8222 1880420.9
190 1441.8052 -3259.0364 -3163.8023 -3196.3556 1904512.1
200 1569.0317 -3265.0089 -3161.3714 -3196.3328 1899462.7
Loop time of 0.348631 on 4 procs for 100 steps with 512 atoms
Loop time of 0.098053 on 1 procs for 100 steps with 512 atoms
Performance: 24.783 ns/day, 0.968 hours/ns, 286.836 timesteps/s
81.0% CPU use with 4 MPI tasks x 1 OpenMP threads
Performance: 88.116 ns/day, 0.272 hours/ns, 1019.857 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 | 0.13281 | 0.15657 | 0.20106 | 6.9 | 44.91
Neigh | 0.00037527 | 0.00039309 | 0.00040412 | 0.0 | 0.11
Comm | 0.12177 | 0.16672 | 0.19154 | 6.8 | 47.82
Output | 0.00019097 | 0.000462 | 0.0012722 | 0.0 | 0.13
Modify | 0.018353 | 0.020198 | 0.02302 | 1.3 | 5.79
Other | | 0.004286 | | | 1.23
Pair | 0.096055 | 0.096055 | 0.096055 | 0.0 | 97.96
Neigh | 0.00079703 | 0.00079703 | 0.00079703 | 0.0 | 0.81
Comm | 0.00034523 | 0.00034523 | 0.00034523 | 0.0 | 0.35
Output | 6.8903e-05 | 6.8903e-05 | 6.8903e-05 | 0.0 | 0.07
Modify | 0.00060797 | 0.00060797 | 0.00060797 | 0.0 | 0.62
Other | | 0.0001793 | | | 0.18
Nlocal: 128.000 ave 132 max 123 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Nghost: 529.500 ave 533 max 524 min
Histogram: 1 0 0 0 0 0 1 1 0 1
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1028.00 ave 1028 max 1028 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 3651.00 ave 3783 max 3494 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 14604.0 ave 14604 max 14604 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 14604
Ave neighs/atom = 28.523438
Neighbor list builds = 1
Dangerous builds = 0
# Test Tersoff/ZBL model for SiC
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.equil
Reading restart file ...
restart file = 29 Sep 2021, LAMMPS = 29 Sep 2021
restoring atom style atomic from restart
orthogonal box = (0.0000000 0.0000000 0.0000000) to (21.724000 21.724000 21.724000)
1 by 1 by 1 MPI processor grid
pair style sw stores no restart info
512 atoms
read_restart CPU = 0.000 seconds
variable fac equal 0.6
change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale ${fac} z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale ${fac} remap
change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap
Changing box ...
orthogonal box = (4.3448000 0.0000000 0.0000000) to (17.379200 21.724000 21.724000)
orthogonal box = (4.3448000 4.3448000 0.0000000) to (17.379200 17.379200 21.724000)
orthogonal box = (4.3448000 4.3448000 4.3448000) to (17.379200 17.379200 17.379200)
pair_style tersoff/zbl
pair_coeff * * SiC.tersoff.zbl C C C C Si Si Si Si
Reading tersoff/zbl potential file SiC.tersoff.zbl with DATE: 2009-04-15
thermo_style custom step temp epair etotal econserve press
thermo 10
fix 1 all nvt temp $t $t 0.1
fix 1 all nvt temp 1800 $t 0.1
fix 1 all nvt temp 1800 1800 0.1
Resetting global fix info from restart file:
fix style: nvt, fix ID: 1
timestep 1.0e-3
neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
run 100
All restart file global fix info was re-assigned
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4
ghost atom cutoff = 4
binsize = 2, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair tersoff/zbl, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.002 | 3.002 | 3.002 Mbytes
Step Temp E_pair TotEng Econserve Press
100 1112.9699 7067.9634 7141.4772 7129.3415 17683957
110 1676.669 7033.1458 7143.893 7128.6921 17837566
120 2450.2667 6982.2491 7144.094 7126.9524 18220027
130 2726.9659 6964.1219 7144.2432 7126.7678 18230324
140 2729.421 6962.7393 7143.0228 7127.2074 18176317
150 2738.5449 6959.1761 7140.0623 7127.6671 18068370
160 2687.2419 6958.1183 7135.6158 7127.8492 18156214
170 2697.7325 6952.1482 7130.3387 7127.7898 17978251
180 2577.9885 6954.5611 7124.8422 7127.5615 18068920
190 2502.6928 6954.4558 7119.7635 7127.67 18049652
200 2517.4866 6947.962 7114.2469 7127.1972 18209451
Loop time of 0.810948 on 1 procs for 100 steps with 512 atoms
Performance: 10.654 ns/day, 2.253 hours/ns, 123.312 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.8082 | 0.8082 | 0.8082 | 0.0 | 99.66
Neigh | 0.001195 | 0.001195 | 0.001195 | 0.0 | 0.15
Comm | 0.00054765 | 0.00054765 | 0.00054765 | 0.0 | 0.07
Output | 0.0001018 | 0.0001018 | 0.0001018 | 0.0 | 0.01
Modify | 0.00062656 | 0.00062656 | 0.00062656 | 0.0 | 0.08
Other | | 0.0002768 | | | 0.03
Nlocal: 512.000 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1542.00 ave 1542 max 1542 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0.00000 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 30142.0 ave 30142 max 30142 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 30142
Ave neighs/atom = 58.871094
Neighbor list builds = 1
Dangerous builds = 0
shell rm restart.equil
Total wall time: 0:00:01

View File

@ -17,7 +17,7 @@ parser = ArgumentParser(prog='Install.py',
# settings
version = "2.7.1"
version = "2.7.4"
mode = "static"
# help message
@ -51,9 +51,14 @@ checksums = { \
'2.6.0' : '204d2edae58d9b10ba3ad460cad64191', \
'2.6.1' : '89a9a450fc6025299fe16af235957163', \
'2.6.3' : 'a9f8028fd74528c2024781ea1fdefeee', \
'2.6.5' : 'b67356f027e5c2747823b0422c3b0ec2', \
'2.6.6' : '6b470dcdce04c221ea42d8500b03c49b', \
'2.7.0' : '95f29dd0c067577f11972ff90dfc7d12', \
'2.7.1' : '4eac6a462ec84dfe0cec96c82421b8e8', \
'2.7.2' : 'cfa0b4dd90a81c25d3302e8d97bfeaea', \
'2.7.3' : 'f00cc82edfefe6bb3df934911dbe32fb', \
'2.7.4' : 'f858e0b6aed173748fc85b6bc8a9dcb3', \
'2.8.0' : '489b23daba70da78cf0506cbc31689c6', \
}
# parse and process arguments

View File

@ -1,4 +1,4 @@
# DATE: 2017-11-28 CONTRIBUTOR: J.H. Los, J.M.H. Kroes CITATION: Los et al. Phys. Rev. B 96, 184108 (2017)
# UNITS: metal DATE: 2017-11-28 CONTRIBUTOR: J.H. Los, J.M.H. Kroes CITATION: Los et al. Phys. Rev. B 96, 184108 (2017)
# B and N mixture, parameterized for ExTeP potential
@ -9,15 +9,15 @@
# other quantities are unitless
# format of a single entry (one or more lines):
#I J K m, gamma*, lambda3, c, d, h, n, gamma, lambda2, B, R, D, lambda1, A
B B B 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498959 2.5211820 2768.7363631 2.0 0.2 2.6857244 3376.3350735
N N N 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928
B B N 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498959 2.5211820 2768.7363631 2.0 0.2 2.6857244 3376.3350735
N N B 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928
B N B 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0
B N N 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0
N B B 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0
N B N 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0
#I J K m, gamma*, lambda3, c, d, h, n, gamma, lambda2, B, R, D, lambda1, A
B B B 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498959 2.5211820 2768.7363631 2.0 0.2 2.6857244 3376.3350735
N N N 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928
B B N 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498959 2.5211820 2768.7363631 2.0 0.2 2.6857244 3376.3350735
N N B 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928
B N B 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2.0 0.2 2.95 3330.0655849887
B N N 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2.0 0.2 2.95 3330.0655849887
N B B 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2.0 0.2 2.95 3330.0655849887
N B N 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2.0 0.2 2.95 3330.0655849887
#
# 1.9925 Bicubic Splines Parameters
#

View File

@ -122,6 +122,9 @@ class lammps(object):
for f in os.listdir(winpath)]):
lib_ext = ".dll"
modpath = winpath
elif any([f.startswith('liblammps') and f.endswith('.so')
for f in os.listdir(modpath)]):
lib_ext = ".so"
else:
import platform
if platform.system() == "Darwin":

View File

@ -1,5 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -19,13 +17,14 @@
#include "pair_lj_cubic.h"
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include <cmath>
#include "pair_lj_cubic_const.h"
@ -34,7 +33,7 @@ using namespace PairLJCubicConstants;
/* ---------------------------------------------------------------------- */
PairLJCubic::PairLJCubic(LAMMPS *lmp) : Pair(lmp) {}
PairLJCubic::PairLJCubic(LAMMPS *_lmp) : Pair(_lmp) {}
/* ---------------------------------------------------------------------- */
@ -60,14 +59,14 @@ PairLJCubic::~PairLJCubic()
void PairLJCubic::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj;
double r,t,rmin;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, r2inv, r6inv, forcelj, factor_lj;
double r, t, rmin;
int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -100,40 +99,39 @@ void PairLJCubic::compute(int eflag, int vflag)
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
if (rsq <= cut_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
} else {
r = sqrt(rsq);
rmin = sigma[itype][jtype]*RT6TWO;
t = (r - cut_inner[itype][jtype])/rmin;
forcelj = epsilon[itype][jtype]*(-DPHIDS + A3*t*t/2.0)*r/rmin;
rmin = sigma[itype][jtype] * RT6TWO;
t = (r - cut_inner[itype][jtype]) / rmin;
forcelj = epsilon[itype][jtype] * (-DPHIDS + A3 * t * t / 2.0) * r / rmin;
}
fpair = factor_lj*forcelj*r2inv;
fpair = factor_lj * forcelj * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
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;
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) {
if (rsq <= cut_inner_sq[itype][jtype])
evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]);
else
evdwl = epsilon[itype][jtype]*(PHIS + DPHIDS*t - A3*t*t*t/6.0);
evdwl = epsilon[itype][jtype] * (PHIS + DPHIDS * t - A3 * t * t * t / 6.0);
evdwl *= factor_lj;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
}
@ -149,33 +147,32 @@ void PairLJCubic::compute(int eflag, int vflag)
void PairLJCubic::allocate()
{
allocated = 1;
int n = atom->ntypes;
const int np1 = atom->ntypes + 1;
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(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
memory->create(cut_inner_sq,n+1,n+1,"pair:cut_inner_sq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(cut, np1, np1, "pair:cut");
memory->create(cut_inner, np1, np1, "pair:cut_inner");
memory->create(cut_inner_sq, np1, np1, "pair:cut_inner_sq");
memory->create(epsilon, np1, np1, "pair:epsilon");
memory->create(sigma, np1, np1, "pair:sigma");
memory->create(lj1, np1, np1, "pair:lj1");
memory->create(lj2, np1, np1, "pair:lj2");
memory->create(lj3, np1, np1, "pair:lj3");
memory->create(lj4, np1, np1, "pair:lj4");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJCubic::settings(int narg, char **/*arg*/)
void PairLJCubic::settings(int narg, char ** /*arg*/)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
if (narg != 0) error->all(FLERR, "Illegal pair_style command");
// NOTE: lj/cubic has no global cutoff. instead the cutoff is
// inferred from the lj parameters. so we must not reset cutoffs here.
@ -187,31 +184,30 @@ void PairLJCubic::settings(int narg, char **/*arg*/)
void PairLJCubic::coeff(int narg, char **arg)
{
if (narg != 4)
error->all(FLERR,"Incorrect args for pair coefficients");
if (narg != 4) 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);
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 epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
double rmin = sigma_one*RT6TWO;
double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
double sigma_one = utils::numeric(FLERR, arg[3], false, lmp);
double rmin = sigma_one * RT6TWO;
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut_inner[i][j] = rmin*SS;
cut[i][j] = rmin*SM;
cut_inner[i][j] = rmin * SS;
cut[i][j] = rmin * SM;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -221,22 +217,22 @@ void PairLJCubic::coeff(int narg, char **arg)
double PairLJCubic::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = epsilon[j][i] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = sigma[j][i] = mix_distance(sigma[i][i],sigma[j][j]);
cut_inner[i][j] = cut_inner[j][i] = mix_distance(cut_inner[i][i],
cut_inner[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], sigma[i][i], sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i], sigma[j][j]);
cut_inner[i][j] = mix_distance(cut_inner[i][i], cut_inner[j][j]);
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
}
cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j];
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
cut_inner_sq[i][j] = cut_inner[i][j] * cut_inner[i][j];
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
cut_inner[j][i] = cut_inner[i][j];
cut_inner_sq[j][i] = cut_inner_sq[i][j];
epsilon[j][i] = epsilon[i][j];
sigma[j][i] = sigma[i][j];
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
@ -253,15 +249,15 @@ void PairLJCubic::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
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);
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut_inner[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
fwrite(&epsilon[i][j], sizeof(double), 1, fp);
fwrite(&sigma[i][j], sizeof(double), 1, fp);
fwrite(&cut_inner[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j], sizeof(double), 1, fp);
}
}
}
@ -275,23 +271,23 @@ void PairLJCubic::read_restart(FILE *fp)
read_restart_settings(fp);
allocate();
int i,j;
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 (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,&epsilon[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_inner[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut_inner[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut_inner[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
@ -302,7 +298,7 @@ void PairLJCubic::read_restart(FILE *fp)
void PairLJCubic::write_restart_settings(FILE *fp)
{
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -312,46 +308,42 @@ void PairLJCubic::write_restart_settings(FILE *fp)
void PairLJCubic::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
}
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
if (me == 0) { utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); }
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */
double PairLJCubic::single(int /*i*/, int /*j*/, int itype, int jtype,
double rsq,
double /*factor_coul*/, double factor_lj,
double &fforce)
double PairLJCubic::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj, double &fforce)
{
double r2inv,r6inv,forcelj,philj;
double r,t;
double r2inv, r6inv, forcelj, philj;
double r, t;
double rmin;
// this is a truncated potential with an implicit cutoff
if (rsq >= cutsq[itype][jtype]) {
fforce=0.0;
fforce = 0.0;
return 0.0;
}
r2inv = 1.0/rsq;
r2inv = 1.0 / rsq;
if (rsq <= cut_inner_sq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
} else {
r = sqrt(rsq);
rmin = sigma[itype][jtype]*RT6TWO;
t = (r - cut_inner[itype][jtype])/rmin;
forcelj = epsilon[itype][jtype]*(-DPHIDS + A3*t*t/2.0)*r/rmin;
rmin = sigma[itype][jtype] * RT6TWO;
t = (r - cut_inner[itype][jtype]) / rmin;
forcelj = epsilon[itype][jtype] * (-DPHIDS + A3 * t * t / 2.0) * r / rmin;
}
fforce = factor_lj*forcelj*r2inv;
fforce = factor_lj * forcelj * r2inv;
if (rsq <= cut_inner_sq[itype][jtype])
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]);
else
philj = epsilon[itype][jtype]*(PHIS + DPHIDS*t - A3*t*t*t/6.0);
philj = epsilon[itype][jtype] * (PHIS + DPHIDS * t - A3 * t * t * t / 6.0);
return factor_lj*philj;
return factor_lj * philj;
}

View File

@ -363,13 +363,12 @@ void AngleCharmmIntel::pack_force_const(ForceConst<flt_t> &fc,
template <class flt_t>
void AngleCharmmIntel::ForceConst<flt_t>::set_ntypes(const int nangletypes,
Memory *memory) {
if (memory != nullptr) _memory = memory;
if (nangletypes != _nangletypes) {
if (_nangletypes > 0)
_memory->destroy(fc);
_memory->destroy(fc);
if (nangletypes > 0)
_memory->create(fc,nangletypes,"anglecharmmintel.fc");
}
_nangletypes = nangletypes;
_memory = memory;
}

View File

@ -60,7 +60,7 @@ class AngleCharmmIntel : public AngleCharmm {
} fc_packed1;
fc_packed1 *fc;
ForceConst() : _nangletypes(0) {}
ForceConst() : fc(nullptr), _nangletypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nangletypes, Memory *memory);

View File

@ -343,13 +343,11 @@ void AngleHarmonicIntel::pack_force_const(ForceConst<flt_t> &fc,
template <class flt_t>
void AngleHarmonicIntel::ForceConst<flt_t>::set_ntypes(const int nangletypes,
Memory *memory) {
if (memory != nullptr) _memory = memory;
if (nangletypes != _nangletypes) {
if (_nangletypes > 0)
_memory->destroy(fc);
_memory->destroy(fc);
if (nangletypes > 0)
_memory->create(fc,nangletypes,"anglecharmmintel.fc");
}
_nangletypes = nangletypes;
_memory = memory;
}

View File

@ -60,7 +60,7 @@ class AngleHarmonicIntel : public AngleHarmonic {
} fc_packed1;
fc_packed1 *fc;
ForceConst() : _nangletypes(0) {}
ForceConst() : fc(nullptr), _nangletypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nangletypes, Memory *memory);

View File

@ -323,13 +323,12 @@ void BondFENEIntel::pack_force_const(ForceConst<flt_t> &fc,
template <class flt_t>
void BondFENEIntel::ForceConst<flt_t>::set_ntypes(const int nbondtypes,
Memory *memory) {
if (memory != nullptr) _memory = memory;
if (nbondtypes != _nbondtypes) {
if (_nbondtypes > 0)
_memory->destroy(fc);
_memory->destroy(fc);
if (nbondtypes > 0)
_memory->create(fc,nbondtypes,"bondfeneintel.fc");
}
_nbondtypes = nbondtypes;
_memory = memory;
}

View File

@ -60,7 +60,7 @@ class BondFENEIntel : public BondFENE {
} fc_packed1;
fc_packed1 *fc;
ForceConst() : _nbondtypes(0) {}
ForceConst() : fc(nullptr), _nbondtypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nbondtypes, Memory *memory);

View File

@ -291,13 +291,12 @@ void BondHarmonicIntel::pack_force_const(ForceConst<flt_t> &fc,
template <class flt_t>
void BondHarmonicIntel::ForceConst<flt_t>::set_ntypes(const int nbondtypes,
Memory *memory) {
if (memory != nullptr) _memory = memory;
if (nbondtypes != _nbondtypes) {
if (_nbondtypes > 0)
_memory->destroy(fc);
_memory->destroy(fc);
if (nbondtypes > 0)
_memory->create(fc,nbondtypes,"bondharmonicintel.fc");
}
_nbondtypes = nbondtypes;
_memory = memory;
}

View File

@ -60,7 +60,7 @@ class BondHarmonicIntel : public BondHarmonic {
} fc_packed1;
fc_packed1 *fc;
ForceConst() : _nbondtypes(0) {}
ForceConst() : fc(nullptr), _nbondtypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nbondtypes, Memory *memory);

View File

@ -261,10 +261,10 @@ void DihedralCharmmIntel::eval(const int vflag,
if (c > (flt_t)1.0) c = (flt_t)1.0;
if (c < (flt_t)-1.0) c = (flt_t)-1.0;
const flt_t tcos_shift = fc.bp[type].cos_shift;
const flt_t tsin_shift = fc.bp[type].sin_shift;
const flt_t tk = fc.bp[type].k;
const int m = fc.bp[type].multiplicity;
const flt_t tcos_shift = fc.fc[type].cos_shift;
const flt_t tsin_shift = fc.fc[type].sin_shift;
const flt_t tk = fc.fc[type].k;
const int m = fc.fc[type].multiplicity;
flt_t p = (flt_t)1.0;
flt_t ddf1, df1;
@ -539,10 +539,10 @@ void DihedralCharmmIntel::eval(const int vflag,
(int *) neighbor->dihedrallist[0];
const flt_t * _noalias const weight = &(fc.weight[0]);
const flt_t * _noalias const x_f = &(x[0].x);
const flt_t * _noalias const cos_shift = &(fc.bp[0].cos_shift);
const flt_t * _noalias const sin_shift = &(fc.bp[0].sin_shift);
const flt_t * _noalias const k = &(fc.bp[0].k);
const int * _noalias const multiplicity = &(fc.bp[0].multiplicity);
const flt_t * _noalias const cos_shift = &(fc.fc[0].cos_shift);
const flt_t * _noalias const sin_shift = &(fc.fc[0].sin_shift);
const flt_t * _noalias const k = &(fc.fc[0].k);
const int * _noalias const multiplicity = &(fc.fc[0].multiplicity);
const flt_t * _noalias const plj1 = &(fc.ljp[0][0].lj1);
const flt_t * _noalias const plj2 = &(fc.ljp[0][0].lj2);
const flt_t * _noalias const plj3 = &(fc.ljp[0][0].lj3);
@ -937,8 +937,8 @@ void DihedralCharmmIntel::pack_force_const(ForceConst<flt_t> &fc,
{
const int tp1 = atom->ntypes + 1;
const int bp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(tp1,bp1,memory);
const int dp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(tp1,dp1,memory);
buffers->set_ntypes(tp1);
if (weightflag) {
@ -952,11 +952,11 @@ void DihedralCharmmIntel::pack_force_const(ForceConst<flt_t> &fc,
}
}
for (int i = 1; i < bp1; i++) {
fc.bp[i].multiplicity = multiplicity[i];
fc.bp[i].cos_shift = cos_shift[i];
fc.bp[i].sin_shift = sin_shift[i];
fc.bp[i].k = k[i];
for (int i = 1; i < dp1; i++) {
fc.fc[i].multiplicity = multiplicity[i];
fc.fc[i].cos_shift = cos_shift[i];
fc.fc[i].sin_shift = sin_shift[i];
fc.fc[i].k = k[i];
fc.weight[i] = weight[i];
}
}
@ -965,27 +965,24 @@ void DihedralCharmmIntel::pack_force_const(ForceConst<flt_t> &fc,
template <class flt_t>
void DihedralCharmmIntel::ForceConst<flt_t>::set_ntypes(const int npairtypes,
const int nbondtypes,
const int ndihderaltypes,
Memory *memory) {
if (memory != nullptr) _memory = memory;
if (npairtypes != _npairtypes) {
if (_npairtypes > 0)
_memory->destroy(ljp);
_memory->destroy(ljp);
if (npairtypes > 0)
memory->create(ljp,npairtypes,npairtypes,"fc.ljp");
_memory->create(ljp,npairtypes,npairtypes,"fc.ljp");
}
if (nbondtypes != _nbondtypes) {
if (_nbondtypes > 0) {
_memory->destroy(bp);
_memory->destroy(weight);
}
if (ndihderaltypes != _ndihderaltypes) {
_memory->destroy(fc);
_memory->destroy(weight);
if (nbondtypes > 0) {
_memory->create(bp,nbondtypes,"dihedralcharmmintel.bp");
_memory->create(weight,nbondtypes,"dihedralcharmmintel.weight");
if (ndihderaltypes > 0) {
_memory->create(fc,ndihderaltypes,"dihedralcharmmintel.fc");
_memory->create(weight,ndihderaltypes,"dihedralcharmmintel.weight");
}
}
_npairtypes = npairtypes;
_nbondtypes = nbondtypes;
_memory = memory;
_ndihderaltypes = ndihderaltypes;
}

View File

@ -64,16 +64,16 @@ class DihedralCharmmIntel : public DihedralCharmm {
} fc_packed3;
fc_packed1 **ljp;
fc_packed3 *bp;
fc_packed3 *fc;
flt_t *weight;
ForceConst() : _npairtypes(0), _nbondtypes(0) {}
ForceConst() : ljp(nullptr), fc(nullptr), _npairtypes(0), _ndihderaltypes(0) {}
~ForceConst() { set_ntypes(0, 0, nullptr); }
void set_ntypes(const int npairtypes, const int nbondtypes, Memory *memory);
void set_ntypes(const int npairtypes, const int ndihderaltypes, Memory *memory);
private:
int _npairtypes, _nbondtypes;
int _npairtypes, _ndihderaltypes;
Memory *_memory;
};
ForceConst<float> force_const_single;

View File

@ -225,10 +225,10 @@ void DihedralFourierIntel::eval(const int vflag,
if (EFLAG) deng = (flt_t)0.0;
for (int j = 0; j < nterms[type]; j++) {
const flt_t tcos_shift = fc.bp[j][type].cos_shift;
const flt_t tsin_shift = fc.bp[j][type].sin_shift;
const flt_t tk = fc.bp[j][type].k;
const int m = fc.bp[j][type].multiplicity;
const flt_t tcos_shift = fc.fc[j][type].cos_shift;
const flt_t tsin_shift = fc.fc[j][type].sin_shift;
const flt_t tk = fc.fc[j][type].k;
const int m = fc.fc[j][type].multiplicity;
flt_t p = (flt_t)1.0;
flt_t ddf1, df1;
@ -394,16 +394,16 @@ template <class flt_t, class acc_t>
void DihedralFourierIntel::pack_force_const(ForceConst<flt_t> &fc,
IntelBuffers<flt_t,acc_t> * /*buffers*/)
{
const int bp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(bp1, setflag, nterms, memory);
const int dp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(dp1, setflag, nterms, memory);
for (int i = 1; i < bp1; i++) {
for (int i = 1; i < dp1; i++) {
if (setflag[i]) {
for (int j = 0; j < nterms[i]; j++) {
fc.bp[j][i].cos_shift = cos_shift[i][j];
fc.bp[j][i].sin_shift = sin_shift[i][j];
fc.bp[j][i].k = k[i][j];
fc.bp[j][i].multiplicity = multiplicity[i][j];
fc.fc[j][i].cos_shift = cos_shift[i][j];
fc.fc[j][i].sin_shift = sin_shift[i][j];
fc.fc[j][i].k = k[i][j];
fc.fc[j][i].multiplicity = multiplicity[i][j];
}
}
}
@ -412,22 +412,20 @@ void DihedralFourierIntel::pack_force_const(ForceConst<flt_t> &fc,
/* ---------------------------------------------------------------------- */
template <class flt_t>
void DihedralFourierIntel::ForceConst<flt_t>::set_ntypes(const int nbondtypes,
void DihedralFourierIntel::ForceConst<flt_t>::set_ntypes(const int ndihedraltypes,
int *setflag,
int *nterms,
Memory *memory) {
if (nbondtypes != _nbondtypes) {
if (_nbondtypes > 0)
_memory->destroy(bp);
if (memory != nullptr) _memory = memory;
if (ndihedraltypes != _ndihedraltypes) {
_memory->destroy(fc);
if (nbondtypes > 0) {
if (ndihedraltypes > 0) {
_maxnterms = 1;
for (int i = 1; i <= nbondtypes; i++)
for (int i = 1; i < ndihedraltypes; i++)
if (setflag[i]) _maxnterms = MAX(_maxnterms, nterms[i]);
_memory->create(bp, _maxnterms, nbondtypes, "dihedralfourierintel.bp");
_memory->create(fc, _maxnterms, ndihedraltypes, "dihedralfourierintel.fc");
}
}
_nbondtypes = nbondtypes;
_memory = memory;
_ndihedraltypes = ndihedraltypes;
}

View File

@ -60,15 +60,15 @@ class DihedralFourierIntel : public DihedralFourier {
int multiplicity;
} fc_packed1;
fc_packed1 **bp;
fc_packed1 **fc;
ForceConst() : _nbondtypes(0) {}
ForceConst() : fc(nullptr), _ndihedraltypes(0) {}
~ForceConst() { set_ntypes(0, nullptr, nullptr, nullptr); }
void set_ntypes(const int nbondtypes, int *setflag, int *nterms, Memory *memory);
void set_ntypes(const int ndihedraltypes, int *setflag, int *nterms, Memory *memory);
private:
int _nbondtypes, _maxnterms;
int _ndihedraltypes, _maxnterms;
Memory *_memory;
};
ForceConst<float> force_const_single;

View File

@ -220,10 +220,10 @@ void DihedralHarmonicIntel::eval(const int vflag,
if (c > (flt_t)1.0) c = (flt_t)1.0;
if (c < (flt_t)-1.0) c = (flt_t)-1.0;
const flt_t tcos_shift = fc.bp[type].cos_shift;
const flt_t tsin_shift = fc.bp[type].sin_shift;
const flt_t tk = fc.bp[type].k;
const int m = fc.bp[type].multiplicity;
const flt_t tcos_shift = fc.fc[type].cos_shift;
const flt_t tsin_shift = fc.fc[type].sin_shift;
const flt_t tk = fc.fc[type].k;
const int m = fc.fc[type].multiplicity;
flt_t p = (flt_t)1.0;
flt_t ddf1, df1;
@ -389,29 +389,28 @@ template <class flt_t, class acc_t>
void DihedralHarmonicIntel::pack_force_const(ForceConst<flt_t> &fc,
IntelBuffers<flt_t,acc_t> * /*buffers*/)
{
const int bp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(bp1,memory);
const int dp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(dp1,memory);
for (int i = 1; i < bp1; i++) {
fc.bp[i].multiplicity = multiplicity[i];
fc.bp[i].cos_shift = cos_shift[i];
fc.bp[i].sin_shift = sin_shift[i];
fc.bp[i].k = k[i];
for (int i = 1; i < dp1; i++) {
fc.fc[i].multiplicity = multiplicity[i];
fc.fc[i].cos_shift = cos_shift[i];
fc.fc[i].sin_shift = sin_shift[i];
fc.fc[i].k = k[i];
}
}
/* ---------------------------------------------------------------------- */
template <class flt_t>
void DihedralHarmonicIntel::ForceConst<flt_t>::set_ntypes(const int nbondtypes,
void DihedralHarmonicIntel::ForceConst<flt_t>::set_ntypes(const int ndihderaltypes,
Memory *memory) {
if (nbondtypes != _nbondtypes) {
if (_nbondtypes > 0)
_memory->destroy(bp);
if (memory != nullptr) _memory = memory;
if (ndihderaltypes != _ndihderaltypes) {
_memory->destroy(fc);
if (nbondtypes > 0)
_memory->create(bp,nbondtypes,"dihedralcharmmintel.bp");
if (ndihderaltypes > 0)
_memory->create(fc,ndihderaltypes,"dihedralcharmmintel.fc");
}
_nbondtypes = nbondtypes;
_memory = memory;
_ndihderaltypes = ndihderaltypes;
}

View File

@ -60,15 +60,15 @@ class DihedralHarmonicIntel : public DihedralHarmonic {
int multiplicity;
} fc_packed1;
fc_packed1 *bp;
fc_packed1 *fc;
ForceConst() : _nbondtypes(0) {}
ForceConst() : fc(nullptr), _ndihderaltypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nbondtypes, Memory *memory);
void set_ntypes(const int ndihderaltypes, Memory *memory);
private:
int _nbondtypes;
int _ndihderaltypes;
Memory *_memory;
};
ForceConst<float> force_const_single;

View File

@ -264,14 +264,14 @@ void DihedralOPLSIntel::eval(const int vflag,
const flt_t sin_4phim = (flt_t)2.0 * cos_2phi * sin_2phim;
flt_t p, pd;
p = fc.bp[type].k1*((flt_t)1.0 + c) +
fc.bp[type].k2*((flt_t)1.0 - cos_2phi) +
fc.bp[type].k3*((flt_t)1.0 + cos_3phi) +
fc.bp[type].k4*((flt_t)1.0 - cos_4phi) ;
pd = fc.bp[type].k1 -
(flt_t)2.0 * fc.bp[type].k2 * sin_2phim +
(flt_t)3.0 * fc.bp[type].k3 * sin_3phim -
(flt_t)4.0 * fc.bp[type].k4 * sin_4phim;
p = fc.fc[type].k1*((flt_t)1.0 + c) +
fc.fc[type].k2*((flt_t)1.0 - cos_2phi) +
fc.fc[type].k3*((flt_t)1.0 + cos_3phi) +
fc.fc[type].k4*((flt_t)1.0 - cos_4phi) ;
pd = fc.fc[type].k1 -
(flt_t)2.0 * fc.fc[type].k2 * sin_2phim +
(flt_t)3.0 * fc.fc[type].k3 * sin_3phim -
(flt_t)4.0 * fc.fc[type].k4 * sin_4phim;
flt_t edihed;
if (EFLAG) edihed = p;
@ -409,29 +409,28 @@ template <class flt_t, class acc_t>
void DihedralOPLSIntel::pack_force_const(ForceConst<flt_t> &fc,
IntelBuffers<flt_t,acc_t> * /*buffers*/)
{
const int bp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(bp1,memory);
const int dp1 = atom->ndihedraltypes + 1;
fc.set_ntypes(dp1,memory);
for (int i = 1; i < bp1; i++) {
fc.bp[i].k1 = k1[i];
fc.bp[i].k2 = k2[i];
fc.bp[i].k3 = k3[i];
fc.bp[i].k4 = k4[i];
for (int i = 1; i < dp1; i++) {
fc.fc[i].k1 = k1[i];
fc.fc[i].k2 = k2[i];
fc.fc[i].k3 = k3[i];
fc.fc[i].k4 = k4[i];
}
}
/* ---------------------------------------------------------------------- */
template <class flt_t>
void DihedralOPLSIntel::ForceConst<flt_t>::set_ntypes(const int nbondtypes,
void DihedralOPLSIntel::ForceConst<flt_t>::set_ntypes(const int ndihderaltypes,
Memory *memory) {
if (nbondtypes != _nbondtypes) {
if (_nbondtypes > 0)
_memory->destroy(bp);
if (memory != nullptr) _memory = memory;
if (ndihderaltypes != _ndihderaltypes) {
_memory->destroy(fc);
if (nbondtypes > 0)
_memory->create(bp,nbondtypes,"dihedralcharmmintel.bp");
if (ndihderaltypes > 0)
_memory->create(fc,ndihderaltypes,"dihedralcharmmintel.fc");
}
_nbondtypes = nbondtypes;
_memory = memory;
_ndihderaltypes = ndihderaltypes;
}

View File

@ -59,15 +59,15 @@ class DihedralOPLSIntel : public DihedralOPLS {
flt_t k1, k2, k3, k4;
} fc_packed1;
fc_packed1 *bp;
fc_packed1 *fc;
ForceConst() : _nbondtypes(0) {}
ForceConst() : fc(nullptr), _ndihderaltypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nbondtypes, Memory *memory);
void set_ntypes(const int ndihderaltypes, Memory *memory);
private:
int _nbondtypes;
int _ndihderaltypes;
Memory *_memory;
};
ForceConst<float> force_const_single;

View File

@ -456,7 +456,7 @@ void FixIntel::pair_init_check(const bool cdmessage)
#endif
int need_tag = 0;
if (atom->molecular != Atom::ATOMIC) need_tag = 1;
if (atom->molecular != Atom::ATOMIC || three_body_neighbor()) need_tag = 1;
// Clear buffers used for pair style
char kmode[80];

View File

@ -426,10 +426,10 @@ template <class flt_t, class acc_t>
void ImproperCvffIntel::pack_force_const(ForceConst<flt_t> &fc,
IntelBuffers<flt_t,acc_t> * /*buffers*/)
{
const int bp1 = atom->nimpropertypes + 1;
fc.set_ntypes(bp1,memory);
const int ip1 = atom->nimpropertypes + 1;
fc.set_ntypes(ip1,memory);
for (int i = 1; i < bp1; i++) {
for (int i = 1; i < ip1; i++) {
fc.fc[i].k = k[i];
fc.fc[i].sign = sign[i];
fc.fc[i].multiplicity = multiplicity[i];
@ -439,15 +439,14 @@ void ImproperCvffIntel::pack_force_const(ForceConst<flt_t> &fc,
/* ---------------------------------------------------------------------- */
template <class flt_t>
void ImproperCvffIntel::ForceConst<flt_t>::set_ntypes(const int nimproper,
void ImproperCvffIntel::ForceConst<flt_t>::set_ntypes(const int nimpropertypes,
Memory *memory) {
if (nimproper != _nimpropertypes) {
if (_nimpropertypes > 0)
_memory->destroy(fc);
if (memory != nullptr) _memory = memory;
if (nimpropertypes != _nimpropertypes) {
_memory->destroy(fc);
if (nimproper > 0)
_memory->create(fc,nimproper,"improperharmonicintel.fc");
if (nimpropertypes > 0)
_memory->create(fc,nimpropertypes,"improperharmonicintel.fc");
}
_nimpropertypes = nimproper;
_memory = memory;
_nimpropertypes = nimpropertypes;
}

View File

@ -62,7 +62,7 @@ class ImproperCvffIntel : public ImproperCvff {
fc_packed1 *fc;
ForceConst() : _nimpropertypes(0) {}
ForceConst() : fc(nullptr), _nimpropertypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nimpropertypes, Memory *memory);

View File

@ -379,10 +379,10 @@ template <class flt_t, class acc_t>
void ImproperHarmonicIntel::pack_force_const(ForceConst<flt_t> &fc,
IntelBuffers<flt_t,acc_t> * /*buffers*/)
{
const int bp1 = atom->nimpropertypes + 1;
fc.set_ntypes(bp1,memory);
const int ip1 = atom->nimpropertypes + 1;
fc.set_ntypes(ip1,memory);
for (int i = 1; i < bp1; i++) {
for (int i = 1; i < ip1; i++) {
fc.fc[i].k = k[i];
fc.fc[i].chi = chi[i];
}
@ -391,15 +391,14 @@ void ImproperHarmonicIntel::pack_force_const(ForceConst<flt_t> &fc,
/* ---------------------------------------------------------------------- */
template <class flt_t>
void ImproperHarmonicIntel::ForceConst<flt_t>::set_ntypes(const int nimproper,
void ImproperHarmonicIntel::ForceConst<flt_t>::set_ntypes(const int nimpropertypes,
Memory *memory) {
if (nimproper != _nimpropertypes) {
if (_nimpropertypes > 0)
_memory->destroy(fc);
if (memory != nullptr) _memory = memory;
if (nimpropertypes != _nimpropertypes) {
_memory->destroy(fc);
if (nimproper > 0)
_memory->create(fc,nimproper,"improperharmonicintel.fc");
if (nimpropertypes > 0)
_memory->create(fc,nimpropertypes,"improperharmonicintel.fc");
}
_nimpropertypes = nimproper;
_memory = memory;
_nimpropertypes = nimpropertypes;
}

View File

@ -61,7 +61,7 @@ class ImproperHarmonicIntel : public ImproperHarmonic {
fc_packed1 *fc;
ForceConst() : _nimpropertypes(0) {}
ForceConst() : fc(nullptr), _nimpropertypes(0) {}
~ForceConst() { set_ntypes(0, nullptr); }
void set_ntypes(const int nimpropertypes, Memory *memory);

View File

@ -207,8 +207,6 @@ void IntelBuffers<flt_t, acc_t>::free_nmax()
template <class flt_t, class acc_t>
void IntelBuffers<flt_t, acc_t>::_grow_nmax(const int offload_end)
{
if (lmp->atom->molecular) _need_tag = 1;
else _need_tag = 0;
#ifdef _LMP_INTEL_OFFLOAD
free_nmax();
int size = lmp->atom->nmax;

View File

@ -55,7 +55,8 @@ NPairSkipIntel::~NPairSkipIntel() {
void NPairSkipIntel::copy_neighbor_info()
{
NPair::copy_neighbor_info();
if (_full_props) delete []_full_props;
// Only need to set _full_props once; npair object deleted for changes
if (_full_props) return;
_full_props = new int[neighbor->nrequest];
for (int i = 0; i < neighbor->nrequest; i++)
_full_props[i] = neighbor->requests[i]->full;

View File

@ -18,10 +18,15 @@
#include "pair_lj_long_coul_long_intel.h"
#include "fix_intel.h"
#include "modify.h"
#include "neigh_request.h"
#include "suffix.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairLJLongCoulLongIntel::PairLJLongCoulLongIntel(LAMMPS *lmp) :
PairLJLongCoulLong(lmp)
{
@ -30,7 +35,22 @@ PairLJLongCoulLongIntel::PairLJLongCoulLongIntel(LAMMPS *lmp) :
cut_respa = nullptr;
}
/* ---------------------------------------------------------------------- */
PairLJLongCoulLongIntel::~PairLJLongCoulLongIntel()
{
}
/* ---------------------------------------------------------------------- */
void PairLJLongCoulLongIntel::init_style()
{
PairLJLongCoulLong::init_style();
neighbor->find_request(this)->intel = 1;
int ifix = modify->find_fix("package_intel");
if (ifix < 0)
error->all(FLERR,"The 'package intel' command is required for /intel styles");
auto fix = static_cast<FixIntel *>(modify->fix[ifix]);
fix->pair_init_check();
}

View File

@ -33,6 +33,7 @@ class PairLJLongCoulLongIntel : public PairLJLongCoulLong {
public:
PairLJLongCoulLongIntel(class LAMMPS *);
virtual ~PairLJLongCoulLongIntel();
void init_style();
};
} // namespace LAMMPS_NS
#endif

View File

@ -88,6 +88,11 @@ void FixNHKokkos<DeviceType>::init()
template<class DeviceType>
void FixNHKokkos<DeviceType>::setup(int /*vflag*/)
{
// tdof needed by compute_temp_target()
t_current = temperature->compute_scalar();
tdof = temperature->dof;
// t_target is needed by NPH and NPT in compute_scalar()
// If no thermostat or using fix nphug,
// t_target must be defined by other means.
@ -737,4 +742,3 @@ template class FixNHKokkos<LMPDeviceType>;
template class FixNHKokkos<LMPHostType>;
#endif
}

View File

@ -142,7 +142,7 @@ void NBinSSAKokkos<DeviceType>::bin_atoms()
k_gbincount.sync<DeviceType>();
ghosts_per_gbin = 0;
NPairSSAKokkosBinIDGhostsFunctor<DeviceType> f(*this);
Kokkos::parallel_reduce(Kokkos::RangePolicy<LMPDeviceType>(nlocal,nall), f, ghosts_per_gbin);
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType>(nlocal,nall), f, ghosts_per_gbin);
}
// actually bin the ghost atoms

View File

@ -263,7 +263,7 @@ struct PairComputeFunctor {
// Loop over neighbors of one atom without coulomb interaction
// This function is called in parallel
KOKKOS_FUNCTION
void compute_item_team(Kokkos::TeamPolicy<>::member_type team,
void compute_item_team(typename Kokkos::TeamPolicy<device_type>::member_type team,
const NeighListKokkos<device_type> &list, const NoCoulTag&) const {
const int inum = team.league_size();
@ -319,7 +319,7 @@ struct PairComputeFunctor {
// Loop over neighbors of one atom with coulomb interaction
// This function is called in parallel
KOKKOS_FUNCTION
void compute_item_team(Kokkos::TeamPolicy<>::member_type team,
void compute_item_team(typename Kokkos::TeamPolicy<device_type>::member_type team,
const NeighListKokkos<device_type> &list, const CoulTag& ) const {
const int inum = team.league_size();
@ -380,7 +380,7 @@ struct PairComputeFunctor {
// Loop over neighbors of one atom without coulomb interaction
// This function is called in parallel
KOKKOS_FUNCTION
EV_FLOAT compute_item_team_ev(Kokkos::TeamPolicy<>::member_type team,
EV_FLOAT compute_item_team_ev(typename Kokkos::TeamPolicy<device_type>::member_type team,
const NeighListKokkos<device_type> &list, const NoCoulTag&) const {
EV_FLOAT ev;
@ -475,7 +475,7 @@ struct PairComputeFunctor {
// Loop over neighbors of one atom with coulomb interaction
// This function is called in parallel
KOKKOS_FUNCTION
EV_FLOAT compute_item_team_ev(Kokkos::TeamPolicy<>::member_type team,
EV_FLOAT compute_item_team_ev(typename Kokkos::TeamPolicy<device_type>::member_type team,
const NeighListKokkos<device_type> &list, const CoulTag& ) const {
EV_FLOAT ev;
@ -684,12 +684,12 @@ struct PairComputeFunctor {
}
KOKKOS_INLINE_FUNCTION
void operator()(const typename Kokkos::TeamPolicy<>::member_type& team) const {
void operator()(const typename Kokkos::TeamPolicy<device_type>::member_type& team) const {
compute_item_team(team,list,typename DoCoul<PairStyle::COUL_FLAG>::type());
}
KOKKOS_INLINE_FUNCTION
void operator()(const typename Kokkos::TeamPolicy<>::member_type& team, value_type &energy_virial) const {
void operator()(const typename Kokkos::TeamPolicy<device_type>::member_type& team, value_type &energy_virial) const {
energy_virial += compute_item_team_ev(team,list,typename DoCoul<PairStyle::COUL_FLAG>::type());
}
};
@ -711,7 +711,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename std::enable_if<!((NE
return ev;
}
template<class FunctorStyle>
template<class DeviceType, class FunctorStyle>
int GetTeamSize(FunctorStyle& KOKKOS_GPU_ARG(functor), int KOKKOS_GPU_ARG(inum),
int KOKKOS_GPU_ARG(reduce_flag), int team_size, int KOKKOS_GPU_ARG(vector_length)) {
@ -719,9 +719,9 @@ int GetTeamSize(FunctorStyle& KOKKOS_GPU_ARG(functor), int KOKKOS_GPU_ARG(inum),
int team_size_max;
if (reduce_flag)
team_size_max = Kokkos::TeamPolicy<>(inum,Kokkos::AUTO).team_size_max(functor,Kokkos::ParallelReduceTag());
team_size_max = Kokkos::TeamPolicy<DeviceType>(inum,Kokkos::AUTO).team_size_max(functor,Kokkos::ParallelReduceTag());
else
team_size_max = Kokkos::TeamPolicy<>(inum,Kokkos::AUTO).team_size_max(functor,Kokkos::ParallelForTag());
team_size_max = Kokkos::TeamPolicy<DeviceType>(inum,Kokkos::AUTO).team_size_max(functor,Kokkos::ParallelForTag());
if (team_size*vector_length > team_size_max)
team_size = team_size_max/vector_length;
@ -746,14 +746,14 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename std::enable_if<(NEIG
if (fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) {
PairComputeFunctor<PairStyle,NEIGHFLAG,false,Specialisation > ff(fpair,list);
atoms_per_team = GetTeamSize(ff, list->inum, (fpair->eflag || fpair->vflag), atoms_per_team, vector_length);
Kokkos::TeamPolicy<Kokkos::IndexType<int> > policy(list->inum,atoms_per_team,vector_length);
atoms_per_team = GetTeamSize<typename PairStyle::device_type>(ff, list->inum, (fpair->eflag || fpair->vflag), atoms_per_team, vector_length);
Kokkos::TeamPolicy<typename PairStyle::device_type,Kokkos::IndexType<int> > policy(list->inum,atoms_per_team,vector_length);
if (fpair->eflag || fpair->vflag) Kokkos::parallel_reduce(policy,ff,ev);
else Kokkos::parallel_for(policy,ff);
} else {
PairComputeFunctor<PairStyle,NEIGHFLAG,true,Specialisation > ff(fpair,list);
atoms_per_team = GetTeamSize(ff, list->inum, (fpair->eflag || fpair->vflag), atoms_per_team, vector_length);
Kokkos::TeamPolicy<Kokkos::IndexType<int> > policy(list->inum,atoms_per_team,vector_length);
atoms_per_team = GetTeamSize<typename PairStyle::device_type>(ff, list->inum, (fpair->eflag || fpair->vflag), atoms_per_team, vector_length);
Kokkos::TeamPolicy<typename PairStyle::device_type,Kokkos::IndexType<int> > policy(list->inum,atoms_per_team,vector_length);
if (fpair->eflag || fpair->vflag) Kokkos::parallel_reduce(policy,ff,ev);
else Kokkos::parallel_for(policy,ff);
}

View File

@ -669,7 +669,7 @@ int PairMultiLucyRXKokkos<DeviceType>::pack_forward_comm_kokkos(int n, DAT::tdua
d_sendlist = k_sendlist.view<DeviceType>();
iswap = iswap_in;
v_buf = buf.view<DeviceType>();
Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagPairMultiLucyRXPackForwardComm>(0,n),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXPackForwardComm>(0,n),*this);
return n;
}
@ -687,7 +687,7 @@ void PairMultiLucyRXKokkos<DeviceType>::unpack_forward_comm_kokkos(int n, int fi
{
first = first_in;
v_buf = buf.view<DeviceType>();
Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagPairMultiLucyRXUnpackForwardComm>(0,n),*this);
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairMultiLucyRXUnpackForwardComm>(0,n),*this);
atomKK->modified(execution_space,DPDRHO_MASK); // needed for auto_sync
}

View File

@ -149,7 +149,6 @@ void PairTersoffKokkos<DeviceType>::setup_params()
}
k_params.template modify<LMPHostType>();
}
/* ---------------------------------------------------------------------- */
@ -171,7 +170,7 @@ void PairTersoffKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.view<DeviceType>();
}
if (vflag_atom) {
if (vflag_either) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.view<DeviceType>();
@ -271,7 +270,7 @@ void PairTersoffKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
if (vflag_either) {
if (need_dup)
Kokkos::Experimental::contribute(d_vatom, dup_vatom);
k_vatom.template modify<DeviceType>();
@ -299,6 +298,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeShortNeigh,
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
const F_FLOAT cutmax_sq = cutmax*cutmax;
const int jnum = d_numneigh[i];
int inside = 0;
@ -311,7 +311,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeShortNeigh,
const X_FLOAT delz = ztmp - x(j,2);
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutmax*cutmax) {
if (rsq < cutmax_sq) {
d_neighbors_short(i,inside) = j;
inside++;
}
@ -480,7 +480,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeHalf<NEIGHFL
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (vflag_atom) {
if (vflag_either) {
F_FLOAT delrij[3], delrik[3];
delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
@ -639,7 +639,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
f_y += fi[1];
f_z += fi[2];
if (vflag_atom) {
if (vflag_either) {
F_FLOAT delrij[3], delrik[3];
delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
@ -764,7 +764,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
f_y += fj[1];
f_z += fj[2];
if (vflag_atom) {
if (vflag_either) {
F_FLOAT delrji[3], delrjk[3];
delrji[0] = -delx1; delrji[1] = -dely1; delrji[2] = -delz1;
delrjk[0] = -delx2; delrjk[1] = -dely2; delrjk[2] = -delz2;
@ -924,10 +924,13 @@ double PairTersoffKokkos<DeviceType>::ters_dbij(const int &i, const int &j,
const int &k, const F_FLOAT &bo) const
{
const F_FLOAT tmp = paramskk(i,j,k).beta * bo;
if (tmp > paramskk(i,j,k).c1) return paramskk(i,j,k).beta * -0.5/sqrt(tmp*tmp);//*pow(tmp,-1.5);
const F_FLOAT factor = -0.5/sqrt(tmp*tmp*tmp); //pow(tmp,-1.5)
if (tmp > paramskk(i,j,k).c1) return paramskk(i,j,k).beta * factor;
if (tmp > paramskk(i,j,k).c2)
return paramskk(i,j,k).beta * (-0.5/sqrt(tmp*tmp) * //*pow(tmp,-1.5) *
(1.0 - 0.5*(1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
return paramskk(i,j,k).beta * (factor *
// error in negligible 2nd term fixed 2/21/2022
// (1.0 - 0.5*(1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
(1.0 - (1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
pow(tmp,-paramskk(i,j,k).powern)));
if (tmp < paramskk(i,j,k).c4) return 0.0;
if (tmp < paramskk(i,j,k).c3)
@ -1127,7 +1130,6 @@ void PairTersoffKokkos<DeviceType>::ters_dthbk(
vec3_scaleadd(fc*dgijk*ex_delr,dcosfk,fk,fk);
vec3_scaleadd(-fc*gijk*dex_delr,rik_hat,fk,fk);
vec3_scale(prefactor,fk,fk);
}
/* ---------------------------------------------------------------------- */
@ -1216,12 +1218,12 @@ void PairTersoffKokkos<DeviceType>::v_tally3(EV_FLOAT &ev, const int &i, const i
F_FLOAT v[6];
v[0] = THIRD * (drij[0]*fj[0] + drik[0]*fk[0]);
v[1] = THIRD * (drij[1]*fj[1] + drik[1]*fk[1]);
v[2] = THIRD * (drij[2]*fj[2] + drik[2]*fk[2]);
v[3] = THIRD * (drij[0]*fj[1] + drik[0]*fk[1]);
v[4] = THIRD * (drij[0]*fj[2] + drik[0]*fk[2]);
v[5] = THIRD * (drij[1]*fj[2] + drik[1]*fk[2]);
v[0] = (drij[0]*fj[0] + drik[0]*fk[0]);
v[1] = (drij[1]*fj[1] + drik[1]*fk[1]);
v[2] = (drij[2]*fj[2] + drik[2]*fk[2]);
v[3] = (drij[0]*fj[1] + drik[0]*fk[1]);
v[4] = (drij[0]*fj[2] + drik[0]*fk[2]);
v[5] = (drij[1]*fj[2] + drik[1]*fk[2]);
if (vflag_global) {
ev.v[0] += v[0];
@ -1233,6 +1235,13 @@ void PairTersoffKokkos<DeviceType>::v_tally3(EV_FLOAT &ev, const int &i, const i
}
if (vflag_atom) {
v[0] *= THIRD;
v[1] *= THIRD;
v[2] *= THIRD;
v[3] *= THIRD;
v[4] *= THIRD;
v[5] *= THIRD;
a_vatom(i,0) += v[0]; a_vatom(i,1) += v[1]; a_vatom(i,2) += v[2];
a_vatom(i,3) += v[3]; a_vatom(i,4) += v[4]; a_vatom(i,5) += v[5];
if (NEIGHFLAG != FULL) {
@ -1242,7 +1251,6 @@ void PairTersoffKokkos<DeviceType>::v_tally3(EV_FLOAT &ev, const int &i, const i
a_vatom(k,3) += v[3]; a_vatom(k,4) += v[4]; a_vatom(k,5) += v[5];
}
}
}
/* ---------------------------------------------------------------------- */

View File

@ -31,6 +31,7 @@
#include "memory_kokkos.h"
#include "error.h"
#include "atom_masks.h"
#include "suffix.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -44,6 +45,7 @@ template<class DeviceType>
PairTersoffZBLKokkos<DeviceType>::PairTersoffZBLKokkos(LAMMPS *lmp) : PairTersoffZBL(lmp)
{
respa_enable = 0;
suffix_flag |= Suffix::KOKKOS;
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
@ -185,7 +187,7 @@ void PairTersoffZBLKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.view<DeviceType>();
}
if (vflag_atom) {
if (vflag_either) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.view<DeviceType>();
@ -285,7 +287,7 @@ void PairTersoffZBLKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
if (vflag_either) {
if (need_dup)
Kokkos::Experimental::contribute(d_vatom, dup_vatom);
k_vatom.template modify<DeviceType>();
@ -524,7 +526,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<N
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (vflag_atom) {
if (vflag_either) {
F_FLOAT delrij[3], delrik[3];
delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
@ -713,7 +715,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
f_y += fi[1];
f_z += fi[2];
if (vflag_atom) {
if (vflag_either) {
F_FLOAT delrij[3], delrik[3];
delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
@ -838,7 +840,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
f_y += fj[1];
f_z += fj[2];
if (vflag_atom) {
if (vflag_either) {
F_FLOAT delrji[3], delrjk[3];
delrji[0] = -delx1; delrji[1] = -dely1; delrji[2] = -delz1;
delrjk[0] = -delx2; delrjk[1] = -dely2; delrjk[2] = -delz2;
@ -1003,7 +1005,9 @@ double PairTersoffZBLKokkos<DeviceType>::ters_dbij(const int &i, const int &j,
if (tmp > paramskk(i,j,k).c1) return paramskk(i,j,k).beta * -0.5*pow(tmp,-1.5);
if (tmp > paramskk(i,j,k).c2)
return paramskk(i,j,k).beta * (-0.5*pow(tmp,-1.5) *
(1.0 - 0.5*(1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
// error in negligible 2nd term fixed 2/21/2022
//(1.0 - 0.5*(1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
(1.0 - (1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
pow(tmp,-paramskk(i,j,k).powern)));
if (tmp < paramskk(i,j,k).c4) return 0.0;
if (tmp < paramskk(i,j,k).c3)
@ -1313,12 +1317,12 @@ void PairTersoffZBLKokkos<DeviceType>::v_tally3(EV_FLOAT &ev, const int &i, cons
F_FLOAT v[6];
v[0] = THIRD * (drij[0]*fj[0] + drik[0]*fk[0]);
v[1] = THIRD * (drij[1]*fj[1] + drik[1]*fk[1]);
v[2] = THIRD * (drij[2]*fj[2] + drik[2]*fk[2]);
v[3] = THIRD * (drij[0]*fj[1] + drik[0]*fk[1]);
v[4] = THIRD * (drij[0]*fj[2] + drik[0]*fk[2]);
v[5] = THIRD * (drij[1]*fj[2] + drik[1]*fk[2]);
v[0] = (drij[0]*fj[0] + drik[0]*fk[0]);
v[1] = (drij[1]*fj[1] + drik[1]*fk[1]);
v[2] = (drij[2]*fj[2] + drik[2]*fk[2]);
v[3] = (drij[0]*fj[1] + drik[0]*fk[1]);
v[4] = (drij[0]*fj[2] + drik[0]*fk[2]);
v[5] = (drij[1]*fj[2] + drik[1]*fk[2]);
if (vflag_global) {
ev.v[0] += v[0];
@ -1330,6 +1334,13 @@ void PairTersoffZBLKokkos<DeviceType>::v_tally3(EV_FLOAT &ev, const int &i, cons
}
if (vflag_atom) {
v[0] *= THIRD;
v[1] *= THIRD;
v[2] *= THIRD;
v[3] *= THIRD;
v[4] *= THIRD;
v[5] *= THIRD;
a_vatom(i,0) += v[0]; a_vatom(i,1) += v[1]; a_vatom(i,2) += v[2];
a_vatom(i,3) += v[3]; a_vatom(i,4) += v[4]; a_vatom(i,5) += v[5];
if (NEIGHFLAG != FULL) {

View File

@ -67,33 +67,7 @@ MSM::MSM(LAMMPS *lmp)
factors[0] = 2;
MPI_Comm_rank(world,&me);
phi1d = dphi1d = nullptr;
nmax = 0;
part2grid = nullptr;
g_direct = nullptr;
g_direct_top = nullptr;
v0_direct = v1_direct = v2_direct = nullptr;
v3_direct = v4_direct = v5_direct = nullptr;
v0_direct_top = v1_direct_top = v2_direct_top = nullptr;
v3_direct_top = v4_direct_top = v5_direct_top = nullptr;
ngrid = nullptr;
alpha = betax = betay = betaz = nullptr;
nx_msm = ny_msm = nz_msm = nullptr;
nxlo_in = nylo_in = nzlo_in = nullptr;
nxhi_in = nyhi_in = nzhi_in = nullptr;
nxlo_out = nylo_out = nzlo_out = nullptr;
nxhi_out = nyhi_out = nzhi_out = nullptr;
delxinv = delyinv = delzinv = nullptr;
qgrid = nullptr;
egrid = nullptr;
v0grid = v1grid = v2grid = v3grid = v4grid = v5grid = nullptr;
peratom_allocate_flag = 0;
scalar_pressure_flag = 1;
@ -311,6 +285,11 @@ double MSM::estimate_total_error()
void MSM::setup()
{
// change_box may trigger MSM::setup() before MSM::init() was called
// error out and request full initialization.
if (!delxinv) error->all(FLERR, "MSM must be fully initialized for this operation");
double *prd;
double a = cutoff;
@ -626,15 +605,19 @@ void MSM::allocate()
gcall->setup(ngcall_buf1,ngcall_buf2);
npergrid = 1;
memory->destroy(gcall_buf1);
memory->destroy(gcall_buf2);
memory->create(gcall_buf1,npergrid*ngcall_buf1,"msm:gcall_buf1");
memory->create(gcall_buf2,npergrid*ngcall_buf2,"msm:gcall_buf2");
// allocate memory for each grid level
for (int n=0; n<levels; n++) {
memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
memory->create3d_offset(qgrid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:qgrid");
memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid");
@ -654,6 +637,8 @@ void MSM::allocate()
gc[n]->setup(ngc_buf1[n],ngc_buf2[n]);
npergrid = 1;
memory->destroy(gc_buf1[n]);
memory->destroy(gc_buf2[n]);
memory->create(gc_buf1[n],npergrid*ngc_buf1[n],"msm:gc_buf1");
memory->create(gc_buf2[n],npergrid*ngc_buf2[n],"msm:gc_buf2");
} else {
@ -835,11 +820,8 @@ void MSM::deallocate_levels()
{
if (world_levels) {
for (int n=0; n < levels; ++n) {
if (qgrid[n])
memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (egrid[n])
memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
if (gc) {
if (gc[n]) {
@ -1106,6 +1088,8 @@ void MSM::set_grid_global()
h_z = 1.0/tmp[2];
}
deallocate_levels();
// find maximum number of levels
levels = MAX(xlevels,ylevels);
@ -1127,7 +1111,6 @@ void MSM::set_grid_global()
if (!domain->nonperiodic) levels -= 1;
deallocate_levels();
allocate_levels();
// find number of grid levels in each direction
@ -1154,33 +1137,19 @@ void MSM::set_grid_global()
error->all(FLERR,"MSM grid is too large");
// compute number of extra grid points needed for non-periodic boundary conditions
// need to always do this, so we can handle the case of switching from periodic
// to non-periodic.
if (domain->nonperiodic) {
alpha[0] = -(order/2 - 1);
betax[0] = nx_msm[0] + (order/2 - 1);
betay[0] = ny_msm[0] + (order/2 - 1);
betaz[0] = nz_msm[0] + (order/2 - 1);
for (int n = 1; n < levels; n++) {
alpha[n] = -((-alpha[n-1]+1)/2) - (order/2 - 1);
betax[n] = ((betax[n-1]+1)/2) + (order/2 - 1);
betay[n] = ((betay[n-1]+1)/2) + (order/2 - 1);
betaz[n] = ((betaz[n-1]+1)/2) + (order/2 - 1);
}
alpha[0] = -(order/2 - 1);
betax[0] = nx_msm[0] + (order/2 - 1);
betay[0] = ny_msm[0] + (order/2 - 1);
betaz[0] = nz_msm[0] + (order/2 - 1);
for (int n = 1; n < levels; n++) {
alpha[n] = -((-alpha[n-1]+1)/2) - (order/2 - 1);
betax[n] = ((betax[n-1]+1)/2) + (order/2 - 1);
betay[n] = ((betay[n-1]+1)/2) + (order/2 - 1);
betaz[n] = ((betaz[n-1]+1)/2) + (order/2 - 1);
}
if (domain->nonperiodic) {
alpha[0] = -(order/2 - 1);
betax[0] = nx_msm[0] + (order/2 - 1);
betay[0] = ny_msm[0] + (order/2 - 1);
betaz[0] = nz_msm[0] + (order/2 - 1);
for (int n = 1; n < levels; n++) {
alpha[n] = -((-alpha[n-1]+1)/2) - (order/2 - 1);
betax[n] = ((betax[n-1]+1)/2) + (order/2 - 1);
betay[n] = ((betay[n-1]+1)/2) + (order/2 - 1);
betaz[n] = ((betaz[n-1]+1)/2) + (order/2 - 1);
}
}
}
/* ----------------------------------------------------------------------

View File

@ -412,6 +412,14 @@ void FixBondCreate::post_integrate()
int *mask = atom->mask;
int *type = atom->type;
if (constrainflag) {
// communicate partner and 1-2 special neighbors
// to correctly handle angle constraints
commflag = 3;
comm->forward_comm_fix(this);
}
neighbor->build_one(list,1);
inum = list->inum;
ilist = list->ilist;

View File

@ -251,6 +251,7 @@ void PairList::settings(int narg, char **arg)
style[npairs] = NONE;
list_parm_t &par = params[npairs];
par.offset = 0.0;
par.id1 = id1;
par.id2 = id2;

View File

@ -29,7 +29,7 @@ action () {
# enforce package dependency
if (test $1 = 1 || test $1 = 2) then
if (test ! -e ../sna.h) then
echo "Must install SNAP package to use ML-IAP package"
echo "Must install ML-SNAP package to use ML-IAP package"
exit 1
fi
fi
@ -69,7 +69,22 @@ elif (test $1 = 0) then
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
elif (test $1 = 2) then
if (test "$(type cythonize 2> /dev/null)" != "" && test -e ../python_impl.cpp) then
if (type cythonize 2>&1 > /dev/null && test -e ../python_impl.cpp) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*-DMLIAP_PYTHON[^ \t]* //g' ../Makefile.package
fi
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
if (test -e ../Makefile.package) then
sed -i -e 's|^PKG_INC =[ \t]*|&-DMLIAP_PYTHON |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/python\/Makefile.mliap_python
' ../Makefile.package.settings
fi
cythonize -3 ../mliap_model_python_couple.pyx
else
rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h

View File

@ -182,7 +182,6 @@ PACKMOST = \
orient \
peri \
plugin \
poems \
qeq \
reaction \
reaxff \
@ -438,7 +437,7 @@ clean:
clean-all:
rm -rf Obj_*
rm style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h
rm -f style_*.h packages_*.h lmpgitversion.h lmpinstalledpkgs.h
clean-%:
@if [ $@ = "clean-serial" ]; \

View File

@ -111,7 +111,6 @@ void PairReaxFFOMP::init_style()
api->system->n = atom->nlocal; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system
api->system->wsize = comm->nprocs;
if (atom->tag_enable == 0)
@ -119,11 +118,6 @@ void PairReaxFFOMP::init_style()
if (force->newton_pair == 0)
error->all(FLERR,"Pair style reaxff/omp requires newton pair on");
// because system->bigN is an int, we cannot have more atoms than MAXSMALLINT
if (atom->natoms > MAXSMALLINT)
error->all(FLERR,"Too many atoms for pair style reaxff/omp");
// need a half neighbor list w/ Newton off and ghost neighbors
// built whenever re-neighboring occurs
@ -153,7 +147,6 @@ void PairReaxFFOMP::setup()
api->system->n = atom->nlocal; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
oldN = api->system->N;
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system
if (api->system->N > nmax) {
memory->destroy(num_nbrs_offset);
@ -234,7 +227,6 @@ void PairReaxFFOMP::compute(int eflag, int vflag)
api->system->n = atom->nlocal; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system
const int nall = api->system->N;
#if defined(_OPENMP)

View File

@ -73,6 +73,8 @@ reaxc_multi_body_omp.cpp
reaxc_nonbonded_omp.cpp
reaxc_torsion_angles_omp.cpp
reaxc_valence_angles_omp.cpp
fix_qeq_reax.cpp
fix_qeq_reax.h
fix_reaxc.cpp
fix_reaxc.h
fix_reaxc_bonds.cpp

View File

@ -251,8 +251,10 @@ FixReaxFFSpecies::~FixReaxFFSpecies()
if (me == 0) fclose(fp);
if (me == 0 && posflag && multipos_opened) fclose(pos);
modify->delete_compute(fmt::format("SPECATOM_{}",id));
modify->delete_fix(fmt::format("SPECBOND_{}",id));
try {
modify->delete_compute(fmt::format("SPECATOM_{}",id));
modify->delete_fix(fmt::format("SPECBOND_{}",id));
} catch (std::exception &) {}
}
/* ---------------------------------------------------------------------- */

View File

@ -90,7 +90,6 @@ PairReaxFF::PairReaxFF(LAMMPS *lmp) : Pair(lmp)
api->system->num_nbrs = 0;
api->system->n = 0; // my atoms
api->system->N = 0; // mine + ghosts
api->system->bigN = 0; // all atoms in the system
api->system->local_cap = 0;
api->system->total_cap = 0;
api->system->my_atoms = nullptr;
@ -343,7 +342,6 @@ void PairReaxFF::init_style()
api->system->n = atom->nlocal; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system
api->system->wsize = comm->nprocs;
if (atom->tag_enable == 0)
@ -351,11 +349,6 @@ void PairReaxFF::init_style()
if (force->newton_pair == 0)
error->all(FLERR,"Pair style reaxff requires newton pair on");
// because system->bigN is an int, we cannot have more atoms than MAXSMALLINT
if (atom->natoms > MAXSMALLINT)
error->all(FLERR,"Too many atoms for pair style reaxff");
// need a half neighbor list w/ Newton off and ghost neighbors
// built whenever re-neighboring occurs
@ -383,7 +376,6 @@ void PairReaxFF::setup()
api->system->n = atom->nlocal; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
oldN = api->system->N;
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system
if (setup_flag == 0) {
@ -464,7 +456,6 @@ void PairReaxFF::compute(int eflag, int vflag)
api->system->n = atom->nlocal; // my atoms
api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
api->system->bigN = static_cast<int> (atom->natoms); // all atoms in the system
// setup data structures

View File

@ -195,7 +195,6 @@ struct LR_lookup_table; // forward declaration
struct reax_system {
reax_interaction reax_param;
rc_bigint bigN;
int n, N, numH;
int local_cap, total_cap, Hcap;
int wsize, my_rank, num_nbrs;

View File

@ -1822,6 +1822,11 @@ void FixRigid::setup_bodies_static()
int *inbody;
if (inpfile) {
// must call it here so it doesn't override read in data but
// initialize bodies whose dynamic settings not set in inpfile
setup_bodies_dynamic();
memory->create(inbody,nbody,"rigid:inbody");
for (ibody = 0; ibody < nbody; ibody++) inbody[ibody] = 0;
readfile(0,masstotal,xcm,vcm,angmom,imagebody,inbody);
@ -2290,11 +2295,17 @@ void FixRigid::readfile(int which, double *vec,
if (*start != '\0' && *start != '#') break;
}
sscanf(line,"%d",&nlines);
if (sscanf(line,"%d",&nlines) != 1) nlines = -1;
if (nlines == 0) fclose(fp);
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
if (nlines == 0) error->all(FLERR,"Fix rigid file has no lines");
// empty file with 0 lines is needed to trigger initial restart file
// generation when no infile was previously used.
if (nlines == 0) return;
else if (nlines < 0) error->all(FLERR,"Fix rigid file has incorrect format");
char *buffer = new char[CHUNK*MAXLINE];
char **values = new char*[ATTRIBUTE_PERBODY];
@ -2406,7 +2417,7 @@ void FixRigid::write_restart_file(const char *file)
int id;
for (int i = 0; i < nbody; i++) {
if (rstyle == SINGLE || rstyle == GROUP) id = i;
if (rstyle == SINGLE || rstyle == GROUP) id = i+1;
else id = body2mol[i];
MathExtra::col2mat(ex_space[i],ey_space[i],ez_space[i],p);

View File

@ -1192,13 +1192,7 @@ void FixRigidNHSmall::compute_dof()
for (int k = 0; k < dimension; k++)
if (fabs(b->inertia[k]) < EPSILON) nf_r--;
}
} else if (dimension == 2) {
nf_r = nlocal_body;
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
if (fabs(b->inertia[2]) < EPSILON) nf_r--;
}
}
} else if (dimension == 2) nf_r = nlocal_body;
double nf[2], nfall[2];
nf[0] = nf_t;

View File

@ -1970,6 +1970,11 @@ void FixRigidSmall::setup_bodies_static()
int *inbody;
if (inpfile) {
// must call it here so it doesn't override read in data but
// initialize bodies whose dynamic settings not set in inpfile
setup_bodies_dynamic();
memory->create(inbody,nlocal_body,"rigid/small:inbody");
for (ibody = 0; ibody < nlocal_body; ibody++) inbody[ibody] = 0;
readfile(0,nullptr,inbody);
@ -2463,11 +2468,18 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody)
if (*start != '\0' && *start != '#') break;
}
sscanf(line,"%d",&nlines);
if (sscanf(line,"%d",&nlines) != 1) nlines = -1;
if (nlines == 0) fclose(fp);
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
// empty file with 0 lines is needed to trigger initial restart file
// generation when no infile was previously used.
if (nlines == 0) return;
else if (nlines < 0) error->all(FLERR,"Fix rigid file has incorrect format");
char *buffer = new char[CHUNK*MAXLINE];
char **values = new char*[ATTRIBUTE_PERBODY];

View File

@ -1201,7 +1201,7 @@ void Input::shell()
} else if (strcmp(arg[0],"mkdir") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell mkdir command");
if (me == 0)
if (me == 0) {
for (int i = 1; i < narg; i++) {
#if defined(_WIN32)
rv = _mkdir(arg[i]);
@ -1214,20 +1214,22 @@ void Input::shell()
delete[] message;
}
}
}
} else if (strcmp(arg[0],"mv") == 0) {
if (narg != 3) error->all(FLERR,"Illegal shell mv command");
rv = (rename(arg[1],arg[2]) < 0) ? errno : 0;
MPI_Reduce(&rv,&err,1,MPI_INT,MPI_MAX,0,world);
if (me == 0 && err != 0) {
char *message = shell_failed_message("mv",err);
error->warning(FLERR,message);
delete[] message;
if (me == 0) {
rv = (rename(arg[1],arg[2]) < 0) ? errno : 0;
if (rv != 0) {
char *message = shell_failed_message("mv",err);
error->warning(FLERR,message);
delete[] message;
}
}
} else if (strcmp(arg[0],"rm") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell rm command");
if (me == 0)
if (me == 0) {
for (int i = 1; i < narg; i++) {
if (unlink(arg[i]) < 0) {
char *message = shell_failed_message("rm",errno);
@ -1235,10 +1237,11 @@ void Input::shell()
delete[] message;
}
}
}
} else if (strcmp(arg[0],"rmdir") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell rmdir command");
if (me == 0)
if (me == 0) {
for (int i = 1; i < narg; i++) {
if (rmdir(arg[i]) < 0) {
char *message = shell_failed_message("rmdir",errno);
@ -1246,6 +1249,7 @@ void Input::shell()
delete[] message;
}
}
}
} else if (strcmp(arg[0],"putenv") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell putenv command");

View File

@ -405,7 +405,7 @@ after calling Py_Finalize().
This function can be called to explicitly clear the Python
environment in case it is safe to do so.
.. versionadded:: TBD
.. versionadded:: 20Sep2021
*See also*
:cpp:func:`lammps_mpi_finalize`, :cpp:func:`lammps_kokkos_finalize`

View File

@ -885,23 +885,28 @@ void ReadData::command(int narg, char **arg)
// restore old styles, when reading with nocoeff flag given
if (coeffflag == 0) {
if (force->pair) delete force->pair;
delete force->pair;
delete[] force->pair_style;
force->pair = saved_pair;
force->pair_style = saved_pair_style;
if (force->bond) delete force->bond;
delete force->bond;
delete[] force->bond_style;
force->bond = saved_bond;
force->bond_style = saved_bond_style;
if (force->angle) delete force->angle;
delete force->angle;
delete[] force->angle_style;
force->angle = saved_angle;
force->angle_style = saved_angle_style;
if (force->dihedral) delete force->dihedral;
delete force->dihedral;
delete[] force->dihedral_style;
force->dihedral = saved_dihedral;
force->dihedral_style = saved_dihedral_style;
if (force->improper) delete force->improper;
delete force->improper;
delete[] force->improper_style;
force->improper = saved_improper;
force->improper_style = saved_improper_style;

View File

@ -1,2 +1,2 @@
#define LAMMPS_VERSION "29 Sep 2021"
#define LAMMPS_UPDATE "Update 2"
#define LAMMPS_UPDATE "Update 3"

View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -30,7 +29,7 @@
using namespace LAMMPS_NS;
enum {REGULAR_MODE, CLASS2_MODE};
enum { REGULAR_MODE, CLASS2_MODE };
/* ----------------------------------------------------------------------
called as write_coeff command in input script
@ -39,121 +38,116 @@ enum {REGULAR_MODE, CLASS2_MODE};
void WriteCoeff::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Write_coeff command before simulation box is defined");
error->all(FLERR, "Write_coeff command before simulation box is defined");
if (narg != 1) error->all(FLERR,"Illegal write_coeff command");
if (narg != 1) error->all(FLERR, "Illegal write_coeff command");
char *file = utils::strdup(fmt::format("tmp.{}",arg[0]));
char *file = utils::strdup(fmt::format("{}.tmp", arg[0]));
// initialize relevant styles
lmp->init();
if (comm->me == 0) {
char str[256], coeff[256];
FILE *one = fopen(file,"wb+");
FILE *one = fopen(file, "wb+");
if (one == nullptr)
error->one(FLERR,"Cannot open coeff file {}: {}",
file, utils::getsyserror());
error->one(FLERR, "Cannot open coeff file {}: {}", file, utils::getsyserror());
if (force->pair && force->pair->writedata) {
fprintf(one,"# pair_style %s\npair_coeff\n",force->pair_style);
fprintf(one, "# pair_style %s\npair_coeff\n", force->pair_style);
force->pair->write_data_all(one);
fprintf(one,"end\n");
fprintf(one, "end\n");
}
if (force->bond && force->bond->writedata) {
fprintf(one,"# bond_style %s\nbond_coeff\n",force->bond_style);
fprintf(one, "# bond_style %s\nbond_coeff\n", force->bond_style);
force->bond->write_data(one);
fprintf(one,"end\n");
fprintf(one, "end\n");
}
if (force->angle && force->angle->writedata) {
fprintf(one,"# angle_style %s\nangle_coeff\n",force->angle_style);
fprintf(one, "# angle_style %s\nangle_coeff\n", force->angle_style);
force->angle->write_data(one);
fprintf(one,"end\n");
fprintf(one, "end\n");
}
if (force->dihedral && force->dihedral->writedata) {
fprintf(one,"# dihedral_style %s\ndihedral_coeff\n",
force->dihedral_style);
fprintf(one, "# dihedral_style %s\ndihedral_coeff\n", force->dihedral_style);
force->dihedral->write_data(one);
fprintf(one,"end\n");
fprintf(one, "end\n");
}
if (force->improper && force->improper->writedata) {
fprintf(one,"# improper_style %s\nimproper_coeff\n",
force->improper_style);
fprintf(one, "# improper_style %s\nimproper_coeff\n", force->improper_style);
force->improper->write_data(one);
fprintf(one,"end\n");
fprintf(one, "end\n");
}
rewind(one);
FILE *two = fopen(file+4,"w");
FILE *two = fopen(arg[0], "w");
if (two == nullptr)
error->one(FLERR,"Cannot open coeff file {}: {}",
file+4, utils::getsyserror());
error->one(FLERR, "Cannot open coeff file {}: {}", arg[0], utils::getsyserror());
fprintf(two,"# LAMMPS coeff file via write_coeff, version %s\n",
lmp->version);
fprintf(two, "# LAMMPS coeff file via write_coeff, version %s\n", lmp->version);
while (1) {
int coeff_mode = REGULAR_MODE;
if (fgets(str,256,one) == nullptr) break;
if (fgets(str, 256, one) == nullptr) break;
// some coeffs need special treatment
if (strstr(str,"class2") != nullptr) {
if (strstr(str,"angle_style") != nullptr)
if (strstr(str, "class2") != nullptr) {
if (strstr(str, "angle_style") != nullptr)
coeff_mode = CLASS2_MODE;
else if (strstr(str,"dihedral_style") != nullptr)
else if (strstr(str, "dihedral_style") != nullptr)
coeff_mode = CLASS2_MODE;
else if (strstr(str,"improper_style") != nullptr)
else if (strstr(str, "improper_style") != nullptr)
coeff_mode = CLASS2_MODE;
}
const char *section = (const char *)"";
fputs(str,two); // style
utils::sfgets(FLERR,str,256,one,file,error); // coeff
const char *section = (const char *) "";
fputs(str, two); // style
utils::sfgets(FLERR, str, 256, one, file, error); // coeff
int n = strlen(str);
strcpy(coeff,str);
coeff[n-1] = '\0';
utils::sfgets(FLERR,str,256,one,file,error);
strcpy(coeff, str);
coeff[n - 1] = '\0';
utils::sfgets(FLERR, str, 256, one, file, error);
while (strcmp(str,"end\n") != 0) {
while (strcmp(str, "end\n") != 0) {
if (coeff_mode == REGULAR_MODE) {
fprintf(two,"%s %s",coeff,str);
utils::sfgets(FLERR,str,256,one,file,error);
fprintf(two, "%s %s", coeff, str);
utils::sfgets(FLERR, str, 256, one, file, error);
} else if (coeff_mode == CLASS2_MODE) {
// class2 angles, dihedrals, and impropers can have
// multiple sections and thus need special treatment
if (strcmp(str,"\n") == 0) {
if (strcmp(str, "\n") == 0) {
// all but the the last section end with an empty line.
// skip it and read and parse the next section title
utils::sfgets(FLERR,str,256,one,file,error);
utils::sfgets(FLERR, str, 256, one, file, error);
if (strcmp(str,"BondBond Coeffs\n") == 0)
section = (const char *)"bb";
else if (strcmp(str,"BondAngle Coeffs\n") ==0)
section = (const char *)"ba";
else if (strcmp(str,"MiddleBondTorsion Coeffs\n") == 0)
section = (const char *)"mbt";
else if (strcmp(str,"EndBondTorsion Coeffs\n") == 0)
section = (const char *)"ebt";
else if (strcmp(str,"AngleTorsion Coeffs\n") == 0)
section = (const char *)"at";
else if (strcmp(str,"AngleAngleTorsion Coeffs\n") == 0)
section = (const char *)"aat";
else if (strcmp(str,"BondBond13 Coeffs\n") == 0)
section = (const char *)"bb13";
else if (strcmp(str,"AngleAngle Coeffs\n") == 0)
section = (const char *)"aa";
if (strcmp(str, "BondBond Coeffs\n") == 0)
section = (const char *) "bb";
else if (strcmp(str, "BondAngle Coeffs\n") == 0)
section = (const char *) "ba";
else if (strcmp(str, "MiddleBondTorsion Coeffs\n") == 0)
section = (const char *) "mbt";
else if (strcmp(str, "EndBondTorsion Coeffs\n") == 0)
section = (const char *) "ebt";
else if (strcmp(str, "AngleTorsion Coeffs\n") == 0)
section = (const char *) "at";
else if (strcmp(str, "AngleAngleTorsion Coeffs\n") == 0)
section = (const char *) "aat";
else if (strcmp(str, "BondBond13 Coeffs\n") == 0)
section = (const char *) "bb13";
else if (strcmp(str, "AngleAngle Coeffs\n") == 0)
section = (const char *) "aa";
// gobble up one more empty line
utils::sfgets(FLERR,str,256,one,file,error);
utils::sfgets(FLERR,str,256,one,file,error);
utils::sfgets(FLERR, str, 256, one, file, error);
utils::sfgets(FLERR, str, 256, one, file, error);
}
// parse type number and skip over it
@ -162,11 +156,11 @@ void WriteCoeff::command(int narg, char **arg)
while ((*p != '\0') && (*p == ' ')) ++p;
while ((*p != '\0') && isdigit(*p)) ++p;
fprintf(two,"%s %d %s %s",coeff,type,section,p);
utils::sfgets(FLERR,str,256,one,file,error);
fprintf(two, "%s %d %s %s", coeff, type, section, p);
utils::sfgets(FLERR, str, 256, one, file, error);
}
}
fputc('\n',two);
fputc('\n', two);
}
fclose(one);
fclose(two);

View File

@ -446,3 +446,31 @@ Zr
1.0
0.85
1.15
Cr
2.493879
1.793835
17.641302
19.60545
8.604593
7.170494
1.551848
1.827556
0.18533
0.277995
-2.022754
0.039608
-0.183611
-2.245972
-2.02
0.00
-0.056517
0.439144
0.456
-2.020038
24
51.9961
0.439144
7.170494
0.277995
0.85
1.15

View File

@ -1,22 +1,52 @@
EAM database tool
Xiaowang Zhou (Sandia), xzhou at sandia.gov
based on this paper:
Fortran version (create.f) by Xiaowang Zhou (Sandia), xzhou at sandia.gov
with revisions by Lucas Hale lucas.hale at nist.gov from https://www.ctcms.nist.gov/potentials/entry/2004--Zhou-X-W-Johnson-R-A-Wadley-H-N-G--Al/
Python version (create_eam.py) by Germain Clavier germain.clavier at gmail.com
Most parameters based on this paper:
X. W. Zhou, R. A. Johnson, and H. N. G. Wadley, Phys. Rev. B, 69,
144113 (2004).
Parameters for Cr were taken from:
Lin Z B, Johnson R A and Zhigilei L V, Phys. Rev. B 77 214108 (2008)
This tool can be used to create an DYNAMO-formatted EAM
setfl file for alloy systems, using any combination
setfl file for alloy systems, using any combination
of the elements discussed in the paper and listed in
the EAM_code file, namely:
Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti, Zr
Cu, Ag, Au, Ni, Pd, Pt, Al, Pb, Fe, Mo, Ta, W, Mg, Co, Ti, Zr, Cr
Steps:
WARNING: Please note that the parameter sets used here are all optimized
for the pure metals of the individual elements and that mixing rules will
be applied for creating the inter-element interactions. Those are inferior
to models where the mixed terms were specifically optimized for particular
alloys. Thus any potential files created with this tool should be used
with care and test calculations (e.g. on multiple binary mixtures) performed
to gauge the error.
Steps (create.f):
1) compile create.f -> a.out (e.g. gfortran create.f)
2) edit the input file EAM.input to list 2 or more desired elements to include
3) a.out < EAM.input will create an EAM *.set file
4) in DYNAMO or LAMMPS lingo, this is a setfl file
that can be used with the LAMMPS pair_style eam/alloy command
3) a.out < EAM.input will create an *.eam.alloy potential file
Steps (create_eam.py):
Usage: create_eam.py [-h] [-n NAME [NAME ...]] [-nr NR] [-nrho NRHO]
options:
-n NAME [NAME ...], --names NAME [NAME ...]
Element names
-nr NR Number of point in r space [default 2000]
-nrho NRHO Number of point in rho space [default 2000]
1) you must have numpy installed
2) run "python create_eam.py -n Ta Cu" with the list of desired elements
3) this will create an *.eam.alloy potential file
in DYNAMO or LAMMPS context the created file is referred to as a setfl file
that can be used with the LAMMPS pair_style eam/alloy command

View File

@ -1,4 +1,5 @@
C author: X. W. Zhou, xzhou@sandia.gov
C updates by: Lucas Hale lucas.hale@nist.gov
c open(unit=5,file='a.i')
call inter
c close(5)
@ -9,6 +10,8 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c main subroutine. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine inter
implicit real*8 (a-h,o-z)
implicit integer (i-m)
character*80 atomtype,atommatch,outfile,outelem
namelist /funccard/ atomtype
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
@ -17,9 +20,9 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
common /pass2/ ielement(16),amass(16),Fr(5000,16),
* rhor(5000,16),z2r(5000,16,16),ntypes,blat(16),
* nrho,drho,nr,dr,rc,outfile,outelem
common /pass2/ amass(16),Fr(5000,16),rhor(5000,16),
* z2r(5000,16,16),blat(16),drho,dr,rc,outfile,outelem
common /pass3/ ielement(16),ntypes,nrho,nr
ntypes=0
10 continue
atomtype='none'
@ -69,7 +72,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
read(10,*) ramda1(ntypes)
read(10,*) rhol(ntypes)
read(10,*) rhoh(ntypes)
blat(ntypes)=sqrt(2.0)*re(ntypes)
blat(ntypes)=sqrt(2.0_8)*re(ntypes)
rhoin(ntypes)=rhol(ntypes)*rhoe(ntypes)
rhoout(ntypes)=rhoh(ntypes)*rhoe(ntypes)
else
@ -91,15 +94,15 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if (alatmax .lt. blat(i)) alatmax=blat(i)
if (rhoemax .lt. rhoe(i)) rhoemax=rhoe(i)
2 continue
rc=sqrt(10.0)/2.0*alatmax
rst=0.5
dr=rc/(nr-1.0)
fmax=-1.0
rc=sqrt(10.0_8)/2.0_8*alatmax
rst=0.5_8
dr=rc/(nr-1.0_8)
fmax=-1.0_8
do i1=1,ntypes
do i2=1,i1
if ( i1 .eq. i2) then
do i=1,nr
r=(i-1.0)*dr
r=(i-1.0_8)*dr
if (r .lt. rst) r=rst
call prof(i1,r,fvalue)
if (fmax .lt. fvalue) fmax=fvalue
@ -109,7 +112,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
end do
else
do i=1,nr
r=(i-1.0)*dr
r=(i-1.0_8)*dr
if (r .lt. rst) r=rst
call pair(i1,i2,r,psi)
z2r(i,i1,i2)=r*psi
@ -119,12 +122,13 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
end do
end do
rhom=fmax
if (rhom .lt. 2.0*rhoemax) rhom=2.0*rhoemax
if (rhom .lt. 100.0) rhom=100.0
drho=rhom/(nrho-1.0)
if (rhom .lt. 2.0_8*rhoemax) rhom=2.0_8*rhoemax
if (rhom .lt. 100.0_8) rhom=100.0_8
drho=rhom/(nrho-1.0_8)
do 6 it=1,ntypes
do 7 i=1,nrho
rhoF=(i-1.0)*drho
rhoF=(i-1)*drho
if (i .eq. 1) rhoF=0.0_8
call embed(it,rhoF,emb)
Fr(i,it)=emb
7 continue
@ -135,20 +139,24 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine calculates the electron density. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine prof(it,r,f)
implicit real*8 (a-h,o-z)
implicit integer (i-m)
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
f=fe(it)*exp(-beta1(it)*(r/re(it)-1.0))
f=f/(1.0+(r/re(it)-ramda1(it))**20)
f=fe(it)*exp(-beta1(it)*(r/re(it)-1.0_8))
f=f/(1.0_8+(r/re(it)-ramda1(it))**20)
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine calculates the pair potential. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine pair(it1,it2,r,psi)
implicit real*8 (a-h,o-z)
implicit integer (i-m)
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
@ -156,25 +164,25 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
if (it1 .eq. it2) then
psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0))
psi1=psi1/(1.0+(r/re(it1)-cai(it1))**20)
psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0))
psi2=psi2/(1.0+(r/re(it1)-ramda(it1))**20)
psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0_8))
psi1=psi1/(1.0_8+(r/re(it1)-cai(it1))**20)
psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0_8))
psi2=psi2/(1.0_8+(r/re(it1)-ramda(it1))**20)
psi=psi1-psi2
else
psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0))
psi1=psi1/(1.0+(r/re(it1)-cai(it1))**20)
psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0))
psi2=psi2/(1.0+(r/re(it1)-ramda(it1))**20)
psi1=A(it1)*exp(-alpha(it1)*(r/re(it1)-1.0_8))
psi1=psi1/(1.0_8+(r/re(it1)-cai(it1))**20)
psi2=B(it1)*exp(-beta(it1)*(r/re(it1)-1.0_8))
psi2=psi2/(1.0_8+(r/re(it1)-ramda(it1))**20)
psia=psi1-psi2
psi1=A(it2)*exp(-alpha(it2)*(r/re(it2)-1.0))
psi1=psi1/(1.0+(r/re(it2)-cai(it2))**20)
psi2=B(it2)*exp(-beta(it2)*(r/re(it2)-1.0))
psi2=psi2/(1.0+(r/re(it2)-ramda(it2))**20)
psi1=A(it2)*exp(-alpha(it2)*(r/re(it2)-1.0_8))
psi1=psi1/(1.0_8+(r/re(it2)-cai(it2))**20)
psi2=B(it2)*exp(-beta(it2)*(r/re(it2)-1.0_8))
psi2=psi2/(1.0_8+(r/re(it2)-ramda(it2))**20)
psib=psi1-psi2
call prof(it1,r,f1)
call prof(it2,r,f2)
psi=0.5*(f2/f1*psia+f1/f2*psib)
psi=0.5_8*(f2/f1*psia+f1/f2*psib)
endif
return
end
@ -182,6 +190,8 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This subroutine calculates the embedding energy. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine embed(it,rho,emb)
implicit real*8 (a-h,o-z)
implicit integer (i-m)
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
* ramda1(16),Fi0(16),Fi1(16),Fi2(16),Fi3(16),
@ -193,18 +203,20 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
else
Fm33=Fm4(it)
endif
if (rho .lt. rhoin(it)) then
if (rho .eq. 0.0_8) then
emb = 0.0_8
else if (rho .lt. rhoin(it)) then
emb=Fi0(it)+
* Fi1(it)*(rho/rhoin(it)-1.0)+
* Fi2(it)*(rho/rhoin(it)-1.0)**2+
* Fi3(it)*(rho/rhoin(it)-1.0)**3
* Fi1(it)*(rho/rhoin(it)-1.0_8)+
* Fi2(it)*(rho/rhoin(it)-1.0_8)**2+
* Fi3(it)*(rho/rhoin(it)-1.0_8)**3
else if (rho .lt. rhoout(it)) then
emb=Fm0(it)+
* Fm1(it)*(rho/rhoe(it)-1.0)+
* Fm2(it)*(rho/rhoe(it)-1.0)**2+
* Fm33*(rho/rhoe(it)-1.0)**3
* Fm1(it)*(rho/rhoe(it)-1.0_8)+
* Fm2(it)*(rho/rhoe(it)-1.0_8)**2+
* Fm33*(rho/rhoe(it)-1.0_8)**3
else
emb=Fn(it)*(1.0-fnn(it)*log(rho/rhos(it)))*
emb=Fn(it)*(1.0_8-fnn(it)*log(rho/rhos(it)))*
* (rho/rhos(it))**fnn(it)
endif
return
@ -213,6 +225,8 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c write out set file. c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine writeset
implicit real*8 (a-h,o-z)
implicit integer (i-m)
character*80 outfile,outelem
common /pass1/ re(16),fe(16),rhoe(16),alpha(16),
* beta(16),beta1(16),A(16),B(16),cai(16),ramda(16),
@ -220,16 +234,21 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
* Fm0(16),Fm1(16),Fm2(16),Fm3(16),Fm4(16),
* fnn(16),Fn(16),rhoin(16),rhoout(16),rhol(16),
* rhoh(16),rhos(16)
common /pass2/ ielement(16),amass(16),Fr(5000,16),
* rhor(5000,16),z2r(5000,16,16),ntypes,blat(16),
* nrho,drho,nr,dr,rc,outfile,outelem
common /pass2/ amass(16),Fr(5000,16),rhor(5000,16),
* z2r(5000,16,16),blat(16),drho,dr,rc,outfile,outelem
common /pass3/ ielement(16),ntypes,nrho,nr
character*80 struc
character(8) date
struc='fcc'
outfile = outfile(1:index(outfile,' ')-1)//'.set'
outfile = outfile(1:index(outfile,' ')-1)//'.eam.alloy'
open(unit=1,file=outfile)
write(1,*)
write(1,*)
write(1,*)
call date_and_time(DATE=date)
write(1,*) 'DATE: ',date(1:4),'-',date(5:6),'-',date(7:8),
* ' UNITS: metal CONTRIBUTOR: Xiaowang Zhou xzhou@sandia.gov and ',
* 'Lucas Hale lucas.hale@nist.gov '
write(1,*) 'CITATION: X. W. Zhou, R. A. Johnson, ',
* 'H. N. G. Wadley, Phys. Rev. B, 69, 144113(2004)'
write(1,*) 'Generated by create.f'
write(1,8)ntypes,outelem
8 format(i5,' ',a24)
write(1,9)nrho,drho,nr,dr,rc

View File

@ -0,0 +1,165 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Python version of the code Zhou04_create_v2.f
original author: X. W. Zhou, xzhou@sandia.gov
based on updates by: Lucas Hale lucas.hale@nist.gov
written by: Germain Clavier g.m.g.c.clavier@tue.nl
This script requires the numpy library.
"""
import sys
import argparse as ap
from datetime import date
import numpy as np
from eamDatabase import Database
def prof(at, r):
atom = Database[at]
f = np.zeros(r.shape)
f = atom.fe * np.exp(-atom.beta1 * (r[r >= 0.5] / atom.re - 1.0))
f = f / (1.0 + (r / atom.re - atom.ramda1) ** 20)
return f
def pair(at1, at2, r):
atom = Database[at1]
psi1 = atom.A * np.exp(-atom.alpha * (r / atom.re - 1.0))
psi1 /= 1.0 + (r / atom.re - atom.cai) ** 20
psi2 = atom.B * np.exp(-atom.beta * (r / atom.re - 1.0))
psi2 /= 1.0 + (r / atom.re - atom.ramda) ** 20
if at1 == at2:
psi = psi1 - psi2
else:
psia = psi1 - psi2
atom2 = Database[at2]
psi1 = atom2.A * np.exp(-atom2.alpha * (r / atom2.re - 1.0))
psi1 /= 1.0 + (r / atom2.re - atom2.cai) ** 20
psi2 = atom2.B * np.exp(-atom2.beta * (r / atom2.re - 1.0))
psi2 /= 1.0 + (r / atom2.re - atom2.ramda) ** 20
psib = psi1 - psi2
prof1 = prof(at1, r)
prof2 = prof(at2, r)
psi = 0.5 * (prof2 / prof1 * psia + prof1 / prof2 * psib)
return psi
def embed(at, rho):
atom = Database[at]
Fm33 = np.zeros(rho.shape)
Fm33[rho < atom.rhoe] = atom.Fm3
Fm33[rho >= atom.rhoe] = atom.Fm4
emb = np.zeros(rho.shape)
for i, r in enumerate(rho):
if r == 0:
emb[i] = 0
elif r < atom.rhoin:
dr = r / atom.rhoin - 1
emb[i] = atom.Fi0 + atom.Fi1 * dr + atom.Fi2 * dr**2 + atom.Fi3 * dr**3
elif r < atom.rhoout:
dr = r / atom.rhoe - 1
emb[i] = atom.Fm0 + atom.Fm1 * dr + atom.Fm2 * dr**2 + Fm33[i] * dr**3
else:
dr = r / atom.rhos
emb[i] = atom.Fn * (1.0 - atom.fnn * np.log(dr)) * dr**atom.fnn
return emb
def write_file(attypes, filename, Fr, rhor, z2r, nrho, drho, nr, dr, rc):
struc = "fcc"
with open(filename, "w") as f:
f.write("DATE: " + date.today().strftime("%Y-%m-%d") + " UNITS: metal")
f.write(" CONTRIBUTOR: Xiaowang Zhou xzhou@sandia.gov, Lucas Hale lucas.hale@nist.gov,")
f.write(" and Germain Clavier g.m.g.c.clavier@tue.nl/germain.clavier@gmail.com\n")
f.write(" CITATION: X. W. Zhou, R. A. Johnson, H. N. G. Wadley, Phys. Rev. B, 69, 144113(2004)\n")
f.write("Generated by create_eam.py\n")
f.write("{:<5d} {:<24}\n".format(len(attypes), " ".join(attypes)))
f.write("{:<5d} {:<24.16e} {:<5d} {:<24.16e} {:<24.16e}\n".format(nrho, drho, nr, dr, rc))
for at in attypes:
atom = Database[at]
f.write(
"{:>5d} {:>15.5f} {:>15.5f} {:>8}\n".format(
atom.ielement, atom.amass, atom.blat, struc
)
)
for i, fr in enumerate(Fr[at]):
f.write(" {:>24.16E}".format(fr))
if not (i + 1) % 5:
f.write("\n")
for i, rho in enumerate(rhor[at]):
f.write(" {:>24.16E}".format(rho))
if not (i + 1) % 5:
f.write("\n")
for n1 in range(len(attypes)):
for n2 in range(n1 + 1):
for i, z in enumerate(z2r[n1, n2]):
f.write(" {:>24.16E}".format(z))
if not (i + 1) % 5:
f.write("\n")
def create_eam(argv=None):
parser = ap.ArgumentParser(description="Script to create EAM alloy potential files.")
parser.add_argument("-n", "--names", dest="name", nargs="+",
help="Element names")
parser.add_argument("-nr", dest="nr", type=int, default=2000,
help="Number of point in r space [default 2000]")
parser.add_argument("-nrho", dest="nrho", type=int, default=2000,
help="Number of point in rho space [default 2000]")
args = parser.parse_args(argv)
if not args.name:
parser.print_help()
sys.exit("")
atnames = args.name
nr = args.nr
nrho = args.nrho
for n in atnames:
try:
Database[n]
except KeyError:
output = "Element {} not found in database.\n".format(n)
valid = "Supported elements are: {}".format(" ".join(Database.keys()))
sys.exit("".join([output, valid]))
ntypes = len(atnames)
outfilename = "".join([*atnames, ".eam.alloy"])
rhor = {}
Fr = {}
alatmax = max([Database[at].blat for at in atnames])
rhoemax = max([Database[at].rhoe for at in atnames])
rc = np.sqrt(10.0) / 2.0 * alatmax
rst = 0.5
r = np.linspace(0.0, rc, num=nr, dtype=np.double)
dr = r[1] - r[0]
r[r < rst] = rst
z2r = np.zeros([ntypes, ntypes, nr])
rhomax = -np.inf
for i, n1 in enumerate(atnames):
for j, n2 in enumerate(atnames):
if j > i:
continue
if i == j:
rhor[n1] = prof(n1, r)
rhomax = max(rhomax,np.max(rhor[n1]))
z2r[i, j, :] = r * pair(n1, n2, r)
else:
z2r[i, j, :] = r * pair(n1, n2, r)
z2r = np.where(z2r, z2r, z2r.transpose((1, 0, 2)))
rhomax = max(rhomax, 2.0 * rhoemax, 100.0)
rho = np.linspace(0.0, rhomax, num=nrho, dtype=np.double)
drho = rho[1] - rho[0]
for i, n1 in enumerate(atnames):
Fr[n1] = embed(n1, rho)
write_file(atnames, outfilename, Fr, rhor, z2r, nrho, drho, nr, dr, rc)
if __name__ == "__main__":
try:
create_eam(sys.argv[1:])
except KeyboardInterrupt as exc:
raise SystemExit("User interruption.") from exc

View File

@ -0,0 +1,640 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Python version of the code Zhou04_create_v2.f
original author: X. W. Zhou, xzhou@sandia.gov
based on updates by: Lucas Hale lucas.hale@nist.gov
written by: Germain Clavier g.m.g.c.clavier@tue.nl
This file contains atom attributes for the EAM model and alloy combination. The
original file is EAM_code. It is designed to be used with create_eam.py script.
To add new contribution, just add new AtType instances.
"""
import math
class AtType:
def __init__(
self,
name,
re,
fe,
rhoe,
rhos,
alpha,
beta,
A,
B,
cai,
ramda,
Fi0,
Fi1,
Fi2,
Fi3,
Fm0,
Fm1,
Fm2,
Fm3,
fnn,
Fn,
ielement,
amass,
Fm4,
beta1,
ramda1,
rhol,
rhoh,
):
self.name = name
self.re = re
self.fe = fe
self.rhoe = rhoe
self.rhos = rhos
self.alpha = alpha
self.beta = beta
self.A = A
self.B = B
self.cai = cai
self.ramda = ramda
self.Fi0 = Fi0
self.Fi1 = Fi1
self.Fi2 = Fi2
self.Fi3 = Fi3
self.Fm0 = Fm0
self.Fm1 = Fm1
self.Fm2 = Fm2
self.Fm3 = Fm3
self.fnn = fnn
self.Fn = Fn
self.ielement = ielement
self.amass = amass
self.Fm4 = Fm4
self.beta1 = beta1
self.ramda1 = ramda1
self.rhol = rhol
self.rhoh = rhoh
self.blat = math.sqrt(2.0) * self.re
self.rhoin = self.rhol * self.rhoe
self.rhoout = self.rhoh * self.rhoe
def __repr__(self):
output = """{}:
re = {}; fe = {}
rhoe = {}; rhos = {};
alpha = {}; beta = {};
A = {}; B = {}
cai = {}; ramda = {}
Fi0 = {}; Fi1 = {}; Fi2 = {}; Fi3 = {}
Fm0 = {}; Fm1 = {}; Fm2 = {}; Fm3 = {}; Fm4 = {}
fnn = {}; Fn = {}
ielement = {}; amass = {}
beta1 = {}; ramda1 = {}
rhol = {}; rhoh = {}
blat = {}; rhoin = {}; rhoout = {}"""
return output.format(
self.name,
self.re,
self.fe,
self.rhoe,
self.rhos,
self.alpha,
self.beta,
self.A,
self.B,
self.cai,
self.ramda,
self.Fi0,
self.Fi1,
self.Fi2,
self.Fi3,
self.Fm0,
self.Fm1,
self.Fm2,
self.Fm3,
self.Fm4,
self.fnn,
self.Fn,
self.ielement,
self.amass,
self.beta1,
self.ramda1,
self.rhol,
self.rhoh,
self.blat,
self.rhoin,
self.rhoout,
)
Database = {}
Database["Cu"] = AtType(
"Cu", # Name
2.556162, # re
1.554485, # fe
21.175871, # rhoe
21.175395, # rhos
8.127620, # alpha
4.334731, # beta
0.396620, # A
0.548085, # B
0.308782, # cai
0.756515, # ramda
-2.170269, # Fi0
-0.263788, # Fi1
1.088878, # Fi2
-0.817603, # Fi3
-2.19, # Fm0
0.00, # Fm1
0.561830, # Fm2
-2.100595, # Fm3
0.310490, # fnn
-2.186568, # Fn
29, # ielement
63.546, # amass
-2.100595, # Fm4
4.334731, # beta1
0.756515, # ramda1
0.85, # rhol
1.15, # rhoh
)
Database["Ag"] = AtType(
"Ag",
2.891814,
1.106232,
14.604100,
14.604144,
9.132010,
4.870405,
0.277758,
0.419611,
0.339710,
0.750758,
-1.729364,
-0.255882,
0.912050,
-0.561432,
-1.75,
0.00,
0.744561,
-1.150650,
0.783924,
-1.748423,
47,
107.8682,
-1.150650,
4.870405,
0.750758,
0.85,
1.15,
)
Database["Au"] = AtType(
"Au",
2.885034,
1.529021,
19.991632,
19.991509,
9.516052,
5.075228,
0.229762,
0.356666,
0.356570,
0.748798,
-2.937772,
-0.500288,
1.601954,
-0.835530,
-2.98,
0.00,
1.706587,
-1.134778,
1.021095,
-2.978815,
79,
196.96654,
-1.134778,
5.075228,
0.748798,
0.85,
1.15,
)
Database["Ni"] = AtType(
"Ni",
2.488746,
2.007018,
27.562015,
27.562031,
8.383453,
4.471175,
0.429046,
0.633531,
0.443599,
0.820658,
-2.693513,
-0.076445,
0.241442,
-2.375626,
-2.70,
0.00,
0.265390,
-0.152856,
0.445470,
-2.7,
28,
58.6934,
-0.152856,
4.471175,
0.820658,
0.85,
1.15,
)
Database["Pd"] = AtType(
"Pd",
2.750897,
1.595417,
21.335246,
21.940073,
8.697397,
4.638612,
0.406763,
0.598880,
0.397263,
0.754799,
-2.321006,
-0.473983,
1.615343,
-0.231681,
-2.36,
0.00,
1.481742,
-1.675615,
1.130000,
-2.352753,
46,
106.42,
-1.675615,
4.638612,
0.754799,
0.85,
1.15,
)
Database["Pt"] = AtType(
"Pt",
2.771916,
2.336509,
33.367564,
35.205357,
7.105782,
3.789750,
0.556398,
0.696037,
0.385255,
0.770510,
-1.455568,
-2.149952,
0.528491,
1.222875,
-4.17,
0.00,
3.010561,
-2.420128,
1.450000,
-4.145597,
78,
195.08,
-2.420128,
3.789750,
0.770510,
0.25,
1.15,
)
Database["Al"] = AtType(
"Al",
2.863924,
1.403115,
20.418205,
23.195740,
6.613165,
3.527021,
0.314873,
0.365551,
0.379846,
0.759692,
-2.807602,
-0.301435,
1.258562,
-1.247604,
-2.83,
0.00,
0.622245,
-2.488244,
0.785902,
-2.824528,
13,
26.981539,
-2.488244,
3.527021,
0.759692,
0.85,
1.15,
)
Database["Pb"] = AtType(
"Pb",
3.499723,
0.647872,
8.450154,
8.450063,
9.121799,
5.212457,
0.161219,
0.236884,
0.250805,
0.764955,
-1.422370,
-0.210107,
0.682886,
-0.529378,
-1.44,
0.00,
0.702726,
-0.538766,
0.935380,
-1.439436,
82,
207.2,
-0.538766,
5.212457,
0.764955,
0.85,
1.15,
)
Database["Fe"] = AtType(
"Fe",
2.481987,
1.885957,
20.041463,
20.041463,
9.818270,
5.236411,
0.392811,
0.646243,
0.170306,
0.340613,
-2.534992,
-0.059605,
0.193065,
-2.282322,
-2.54,
0.00,
0.200269,
-0.148770,
0.391750,
-2.539945,
26,
55.847,
-0.148770,
5.236411,
0.340613,
0.85,
1.15,
)
Database["Mo"] = AtType(
"Mo",
2.728100,
2.723710,
29.354065,
29.354065,
8.393531,
4.476550,
0.708787,
1.120373,
0.137640,
0.275280,
-3.692913,
-0.178812,
0.380450,
-3.133650,
-3.71,
0.00,
0.875874,
0.776222,
0.790879,
-3.712093,
42,
95.94,
0.776222,
4.476550,
0.275280,
0.85,
1.15,
)
Database["Ta"] = AtType(
"Ta",
2.860082,
3.086341,
33.787168,
33.787168,
8.489528,
4.527748,
0.611679,
1.032101,
0.176977,
0.353954,
-5.103845,
-0.405524,
1.112997,
-3.585325,
-5.14,
0.00,
1.640098,
0.221375,
0.848843,
-5.141526,
73,
180.9479,
0.221375,
4.527748,
0.353954,
0.85,
1.15,
)
Database["W"] = AtType(
"W",
2.740840,
3.487340,
37.234847,
37.234847,
8.900114,
4.746728,
0.882435,
1.394592,
0.139209,
0.278417,
-4.946281,
-0.148818,
0.365057,
-4.432406,
-4.96,
0.00,
0.661935,
0.348147,
0.582714,
-4.961306,
74,
183.84,
0.348147,
4.746728,
0.278417,
0.85,
1.15,
)
Database["Mg"] = AtType(
"Mg",
3.196291,
0.544323,
7.132600,
7.132600,
10.228708,
5.455311,
0.137518,
0.225930,
0.5,
1.0,
-0.896473,
-0.044291,
0.162232,
-0.689950,
-0.90,
0.00,
0.122838,
-0.226010,
0.431425,
-0.899702,
12,
24.305,
-0.226010,
5.455311,
1.0,
0.85,
1.15,
)
Database["Co"] = AtType(
"Co",
2.505979,
1.975299,
27.206789,
27.206789,
8.679625,
4.629134,
0.421378,
0.640107,
0.5,
1.0,
-2.541799,
-0.219415,
0.733381,
-1.589003,
-2.56,
0.00,
0.705845,
-0.687140,
0.694608,
-2.559307,
27,
58.9332,
-0.687140,
4.629134,
1.0,
0.85,
1.15,
)
Database["Ti"] = AtType(
"Ti",
2.933872,
1.863200,
25.565138,
25.565138,
8.775431,
4.680230,
0.373601,
0.570968,
0.5,
1.0,
-3.203773,
-0.198262,
0.683779,
-2.321732,
-3.22,
0.00,
0.608587,
-0.750710,
0.558572,
-3.219176,
22,
47.88,
-0.750710,
4.680230,
1.0,
0.85,
1.15,
)
Database["Zr"] = AtType(
"Zr",
3.199978,
2.230909,
30.879991,
30.879991,
8.559190,
4.564902,
0.424667,
0.640054,
0.5,
1.0,
-4.485793,
-0.293129,
0.990148,
-3.202516,
-4.51,
0.00,
0.928602,
-0.981870,
0.597133,
-4.509025,
40,
91.224,
-0.981870,
4.564902,
1.0,
0.85,
1.15,
)
Database["Cr"] = AtType(
"Cr",
2.493879,
1.793835,
17.641302,
19.60545,
8.604593,
7.170494,
1.551848,
1.827556,
0.18533,
0.277995,
-2.022754,
0.039608,
-0.183611,
-2.245972,
-2.02,
0.00,
-0.056517,
0.439144,
0.456,
-2.020038,
24,
51.9961,
0.439144,
7.170494,
0.277995,
0.85,
1.15
)

View File

@ -6,14 +6,14 @@ functionality that could be added, built on the Pizza.py modules (as
explained below), send email to Steve Plimpton (sjplimp at
sandia.gov).
log2txt.py convert thermo info in a LAMMPS log file to columns of numbers
logplot.py plot 2 columns of thermo info from a LAMMPS log file
dumpsort.py sort the snapshots of atoms in a LAMMPS dump file by atom ID
dump2cfg.py convert a native LAMMPS dump file to CFG format
dump2xyz.py convert a native LAMMPS dump file to XYZ format
dump2pdb.py convert a native LAMMPS dump file to PDB format
neb_combine.py combine multiple NEB dump files into one time series
neb_final.py combine multiple NEB final states into one sequence of states
log2txt.py convert thermo info in a LAMMPS log file to columns of numbers
logplot.py plot 2 columns of thermo info from a LAMMPS log file
dumpsort.py sort the snapshots of atoms in a LAMMPS dump file by atom ID
dump2cfg.py convert a native LAMMPS dump file to CFG format
dump2xyz.py convert a native LAMMPS dump file to XYZ format
dump2pdb.py convert a native LAMMPS dump file to PDB format
neb_combine.py combine multiple NEB dump files into one time series
neb_final.py combine multiple NEB final states into one sequence of states
See the top of each script file for syntax, or just run it with no
arguments to get a syntax message.
@ -68,7 +68,7 @@ directory where your data files are.
The latter requires 2 things:
1) that the script be made "executable", e.g. type "chmod +x log2txt.py"
2) that the 1st line of the script is the path of the Python installed
2) that the 1st line of the script is the path of the Python installed
on your box, e.g. /usr/local/bin/python
IMPORTANT NOTE: If you run the logplot.py script using the 1st method

View File

@ -16,7 +16,7 @@ from dump import dump
from cfg import cfg
if len(sys.argv) != 8:
raise StandardError, "Syntax: dump2cfg.py dumpfile Nid Ntype Nx Ny Nz cfgfile"
sys.exit("Syntax: dump2cfg.py dumpfile Nid Ntype Nx Ny Nz cfgfile")
dumpfile = sys.argv[1]
nid = int(sys.argv[2])

View File

@ -18,7 +18,7 @@ from dump import dump
from pdbfile import pdbfile
if len(sys.argv) != 8 and len(sys.argv) != 9:
raise StandardError, "Syntax: dump2pdb.py dumpfile Nid Ntype Nx Ny Nz pdbfile template"
sys.exit("Syntax: dump2pdb.py dumpfile Nid Ntype Nx Ny Nz pdbfile template")
dumpfile = sys.argv[1]
nid = int(sys.argv[2])

View File

@ -16,7 +16,7 @@ from dump import dump
from xyz import xyz
if len(sys.argv) != 8:
raise StandardError, "Syntax: dump2xyz.py dumpfile Nid Ntype Nx Ny Nz xyzfile"
sys.exit("Syntax: dump2xyz.py dumpfile Nid Ntype Nx Ny Nz xyzfile")
dumpfile = sys.argv[1]
nid = int(sys.argv[2])

View File

@ -14,7 +14,7 @@ sys.path.append(path)
from dump import dump
if len(sys.argv) != 4:
raise StandardError, "Syntax: dumpsort.py oldfile N newfile"
sys.exit("Syntax: dumpsort.py oldfile N newfile")
oldfile = sys.argv[1]
ncolumn = int(sys.argv[2])

View File

@ -10,23 +10,31 @@
# if no columns listed, all columns are included
# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov
import sys,os
from __future__ import print_function
import sys,os,argparse
path = os.environ["LAMMPS_PYTHON_TOOLS"]
sys.path.append(path)
from log import log
if len(sys.argv) < 3:
raise StandardError, "Syntax: log2txt.py log.lammps data.txt X Y ..."
# set up arg parser
parser = argparse.ArgumentParser()
parser.add_argument('lammpslog', help='name of the lammps log file')
parser.add_argument('outname', help='name of the file to be written')
parser.add_argument('cols', nargs='*', help='any number of column names, optional')
parser.add_argument('-n', action='store_true', help='save column names as the header of the file')
logfile = sys.argv[1]
datafile = sys.argv[2]
columns = sys.argv[3:]
args = parser.parse_args()
logfile = args.lammpslog
datafile = args.outname
columns = args.cols
writenames = args.n
lg = log(logfile)
if columns == []:
lg.write(datafile)
lg.write(datafile, writenames)
else:
str = "lg.write(datafile,"
for word in columns: str += '"' + word + '",'
str = str[:-1] + ')'
str = "lg.write(datafile, %r" % writenames
for word in columns: str += ',"' + word + '"'
str += ')'
eval(str)

View File

@ -1,4 +1,4 @@
#!/usr/bin/env python -i
#!/usr/bin/env python
# Script: logplot.py
# Purpose: use GnuPlot to plot two columns from a LAMMPS log file
@ -8,14 +8,14 @@
# once plot appears, you are in Python interpreter, type C-D to exit
# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov
import sys,os
import sys,os,code
path = os.environ["LAMMPS_PYTHON_TOOLS"]
sys.path.append(path)
from log import log
from gnu import gnu
if len(sys.argv) != 4:
raise StandardError, "Syntax: logplot.py log.lammps X Y"
raise Exception("Syntax: logplot.py log.lammps X Y")
logfile = sys.argv[1]
xlabel = sys.argv[2]
@ -25,4 +25,4 @@ lg = log(logfile)
x,y = lg.get(xlabel,ylabel)
g = gnu()
g.plot(x,y)
print "Type Ctrl-D to exit Python"
code.interact()

View File

@ -38,8 +38,7 @@ while iarg < narg:
else: break
if iarg < narg or not outfile or not rfiles:
print "Syntax: neb_combine.py -o outfile -b backfile -r dump1 dump2 ..."
sys.exit()
sys.exit("Syntax: neb_combine.py -o outfile -b backfile -r dump1 dump2 ...")
if os.path.exists(outfile): os.remove(outfile)

View File

@ -38,8 +38,7 @@ while iarg < narg:
else: break
if iarg < narg or not outfile or not rfiles:
print "Syntax: neb_final.py -o outfile -b backfile -r dump1 dump2 ..."
sys.exit()
sys.exit("Syntax: neb_final.py -o outfile -b backfile -r dump1 dump2 ...")
if os.path.exists(outfile): os.remove(outfile)

View File

@ -6,6 +6,8 @@
# certain rights in this software. This software is distributed under
# the GNU General Public License.
from __future__ import print_function
# cfg tool
oneline = "Convert LAMMPS snapshots to AtomEye CFG format"
@ -65,33 +67,33 @@ class cfg:
ylen = box[4]-box[1]
zlen = box[5]-box[2]
print >>f,"Number of particles = %d " % len(atoms)
print >>f,"# Timestep %d \n#\nA = 1.0 Angstrom" % time
print >>f,"H0(1,1) = %20.10f A " % xlen
print >>f,"H0(1,2) = 0.0 A "
print >>f,"H0(1,3) = 0.0 A "
print >>f,"H0(2,1) = 0.0 A "
print >>f,"H0(2,2) = %20.10f A " % ylen
print >>f,"H0(2,3) = 0.0 A "
print >>f,"H0(3,1) = 0.0 A "
print >>f,"H0(3,2) = 0.0 A "
print >>f,"H0(3,3) = %20.10f A " % zlen
print >>f,"#"
print("Number of particles = %d " % len(atoms), file=f)
print("# Timestep %d \n#\nA = 1.0 Angstrom" % time, file=f)
print("H0(1,1) = %20.10f A " % xlen, file=f)
print("H0(1,2) = 0.0 A ", file=f)
print("H0(1,3) = 0.0 A ", file=f)
print("H0(2,1) = 0.0 A ", file=f)
print("H0(2,2) = %20.10f A " % ylen, file=f)
print("H0(2,3) = 0.0 A ", file=f)
print("H0(3,1) = 0.0 A ", file=f)
print("H0(3,2) = 0.0 A ", file=f)
print("H0(3,3) = %20.10f A " % zlen, file=f)
print("#", file=f)
for atom in atoms:
itype = int(atom[1])
xfrac = (atom[2]-box[0])/xlen
yfrac = (atom[3]-box[1])/ylen
zfrac = (atom[4]-box[2])/zlen
# print >>f,"1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7])
print >>f,"1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac)
# print("1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7]), file=f)
print("1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac), file=f)
print time,
print(time)
sys.stdout.flush()
n += 1
f.close()
print "\nwrote %d snapshots to %s in CFG format" % (n,file)
print("\nwrote %d snapshots to %s in CFG format" % (n,file))
# --------------------------------------------------------------------
@ -120,33 +122,33 @@ class cfg:
ylen = box[4]-box[1]
zlen = box[5]-box[2]
print >>f,"Number of particles = %d " % len(atoms)
print >>f,"# Timestep %d \n#\nA = 1.0 Angstrom" % time
print >>f,"H0(1,1) = %20.10f A " % xlen
print >>f,"H0(1,2) = 0.0 A "
print >>f,"H0(1,3) = 0.0 A "
print >>f,"H0(2,1) = 0.0 A "
print >>f,"H0(2,2) = %20.10f A " % ylen
print >>f,"H0(2,3) = 0.0 A "
print >>f,"H0(3,1) = 0.0 A "
print >>f,"H0(3,2) = 0.0 A "
print >>f,"H0(3,3) = %20.10f A " % zlen
print >>f,"#"
print("Number of particles = %d " % len(atoms), file=f)
print("# Timestep %d \n#\nA = 1.0 Angstrom" % time, file=f)
print("H0(1,1) = %20.10f A " % xlen, file=f)
print("H0(1,2) = 0.0 A ", file=f)
print("H0(1,3) = 0.0 A ", file=f)
print("H0(2,1) = 0.0 A ", file=f)
print("H0(2,2) = %20.10f A " % ylen, file=f)
print("H0(2,3) = 0.0 A ", file=f)
print("H0(3,1) = 0.0 A ", file=f)
print("H0(3,2) = 0.0 A ", file=f)
print("H0(3,3) = %20.10f A " % zlen, file=f)
print("#", file=f)
for atom in atoms:
itype = int(atom[1])
xfrac = (atom[2]-box[0])/xlen
yfrac = (atom[3]-box[1])/ylen
zfrac = (atom[4]-box[2])/zlen
# print >>f,"1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7])
print >>f,"1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac)
# print("1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7]), file=f)
print("1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac), file=f)
print time,
print(time)
sys.stdout.flush()
f.close()
n += 1
print "\nwrote %s snapshots in CFG format" % n
print("\nwrote %s snapshots in CFG format" % n)
# --------------------------------------------------------------------
@ -163,25 +165,25 @@ class cfg:
ylen = box[4]-box[1]
zlen = box[5]-box[2]
print >>f,"Number of particles = %d " % len(atoms)
print >>f,"# Timestep %d \n#\nA = 1.0 Angstrom" % time
print >>f,"H0(1,1) = %20.10f A " % xlen
print >>f,"H0(1,2) = 0.0 A "
print >>f,"H0(1,3) = 0.0 A "
print >>f,"H0(2,1) = 0.0 A "
print >>f,"H0(2,2) = %20.10f A " % ylen
print >>f,"H0(2,3) = 0.0 A "
print >>f,"H0(3,1) = 0.0 A "
print >>f,"H0(3,2) = 0.0 A "
print >>f,"H0(3,3) = %20.10f A " % zlen
print >>f,"#"
print("Number of particles = %d " % len(atoms), file=f)
print("# Timestep %d \n#\nA = 1.0 Angstrom" % time, file=f)
print("H0(1,1) = %20.10f A " % xlen, file=f)
print("H0(1,2) = 0.0 A ", file=f)
print("H0(1,3) = 0.0 A ", file=f)
print("H0(2,1) = 0.0 A ", file=f)
print("H0(2,2) = %20.10f A " % ylen, file=f)
print("H0(2,3) = 0.0 A ", file=f)
print("H0(3,1) = 0.0 A ", file=f)
print("H0(3,2) = 0.0 A ", file=f)
print("H0(3,3) = %20.10f A " % zlen, file=f)
print("#", file=f)
for atom in atoms:
itype = int(atom[1])
xfrac = (atom[2]-box[0])/xlen
yfrac = (atom[3]-box[1])/ylen
zfrac = (atom[4]-box[2])/zlen
# print >>f,"1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7])
print >>f,"1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac)
# print("1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7]), file=f)
print("1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac), file=f)
f.close()

View File

@ -16,15 +16,15 @@ oneline = "Read, write, manipulate dump files and particle attributes"
docstr = """
d = dump("dump.one") read in one or more dump files
d = dump("dump.1 dump.2.gz") can be gzipped
d = dump("dump.*") wildcard expands to multiple files
d = dump("dump.*",0) two args = store filenames, but don't read
d = dump("dump.1 dump.2.gz") can be gzipped
d = dump("dump.*") wildcard expands to multiple files
d = dump("dump.*",0) two args = store filenames, but don't read
incomplete and duplicate snapshots are deleted
if atoms have 5 or 8 columns, assign id,type,x,y,z (ix,iy,iz)
atoms will be unscaled if stored in files as scaled
time = d.next() read next snapshot from dump files
time = d.next() read next snapshot from dump files
used with 2-argument constructor to allow reading snapshots one-at-a-time
snapshot will be skipped only if another snapshot has same time stamp
@ -36,20 +36,20 @@ d.map(1,"id",3,"x") assign names to atom columns (1-N)
not needed if dump file is self-describing
d.tselect.all() select all timesteps
d.tselect.one(N) select only timestep N
d.tselect.none() deselect all timesteps
d.tselect.skip(M) select every Mth step
d.tselect.all() select all timesteps
d.tselect.one(N) select only timestep N
d.tselect.none() deselect all timesteps
d.tselect.skip(M) select every Mth step
d.tselect.test("$t >= 100 and $t < 10000") select matching timesteps
d.delete() delete non-selected timesteps
d.delete() delete non-selected timesteps
selecting a timestep also selects all atoms in the timestep
skip() and test() only select from currently selected timesteps
test() uses a Python Boolean expression with $t for timestep value
Python comparison syntax: == != < > <= >= and or
d.aselect.all() select all atoms in all steps
d.aselect.all(N) select all atoms in one step
d.aselect.all() select all atoms in all steps
d.aselect.all(N) select all atoms in one step
d.aselect.test("$id > 100 and $type == 2") select match atoms in all steps
d.aselect.test("$id > 100 and $type == 2",N) select matching atoms in one step
@ -60,24 +60,24 @@ d.aselect.test("$id > 100 and $type == 2",N) select matching atoms in one step
Python comparison syntax: == != < > <= >= and or
$name must end with a space
d.write("file") write selected steps/atoms to dump file
d.write("file",head,app) write selected steps/atoms to dump file
d.scatter("tmp") write selected steps/atoms to multiple files
d.write("file") write selected steps/atoms to dump file
d.write("file",head,app) write selected steps/atoms to dump file
d.scatter("tmp") write selected steps/atoms to multiple files
write() can be specified with 2 additional flags
headd = 0/1 for no/yes snapshot header, app = 0/1 for write vs append
scatter() files are given timestep suffix: e.g. tmp.0, tmp.100, etc
d.scale() scale x,y,z to 0-1 for all timesteps
d.scale(100) scale atom coords for timestep N
d.unscale() unscale x,y,z to box size to all timesteps
d.unscale(1000) unscale atom coords for timestep N
d.wrap() wrap x,y,z into periodic box via ix,iy,iz
d.unwrap() unwrap x,y,z out of box via ix,iy,iz
d.owrap("other") wrap x,y,z to same image as another atom
d.sort() sort atoms by atom ID in all selected steps
d.sort("x") sort atoms by column value in all steps
d.sort(1000) sort atoms in timestep N
d.scale() scale x,y,z to 0-1 for all timesteps
d.scale(100) scale atom coords for timestep N
d.unscale() unscale x,y,z to box size to all timesteps
d.unscale(1000) unscale atom coords for timestep N
d.wrap() wrap x,y,z into periodic box via ix,iy,iz
d.unwrap() unwrap x,y,z out of box via ix,iy,iz
d.owrap("other") wrap x,y,z to same image as another atom
d.sort() sort atoms by atom ID in all selected steps
d.sort("x") sort atoms by column value in all steps
d.sort(1000) sort atoms in timestep N
scale(), unscale(), wrap(), unwrap(), owrap() operate on all steps and atoms
wrap(), unwrap(), owrap() require ix,iy,iz be defined
@ -89,8 +89,8 @@ d.sort(1000) sort atoms in timestep N
m1,m2 = d.minmax("type") find min/max values for a column
d.set("$ke = $vx * $vx + $vy * $vy") set a column to a computed value
d.setv("type",vector) set a column to a vector of values
d.spread("ke",N,"color") 2nd col = N ints spread over 1st col
d.clone(1000,"color") clone timestep N values to other steps
d.spread("ke",N,"color") 2nd col = N ints spread over 1st col
d.clone(1000,"color") clone timestep N values to other steps
minmax() operates on selected timesteps and atoms
set() operates on selected timesteps and atoms
@ -111,7 +111,7 @@ d.clone(1000,"color") clone timestep N values to other steps
values at every timestep are set to value at timestep N for that atom ID
useful for propagating a color map
t = d.time() return vector of selected timestep values
t = d.time() return vector of selected timestep values
fx,fy,... = d.atom(100,"fx","fy",...) return vector(s) for atom ID N
fx,fy,... = d.vecs(1000,"fx","fy",...) return vector(s) for timestep N
@ -121,8 +121,8 @@ fx,fy,... = d.vecs(1000,"fx","fy",...) return vector(s) for timestep N
index,time,flag = d.iterator(0/1) loop over dump snapshots
time,box,atoms,bonds,tris = d.viz(index) return list of viz objects
d.atype = "color" set column returned as "type" by viz
d.extra("dump.bond") read bond list from dump file
d.extra(data) extract bond/tri/line list from data
d.extra("dump.bond") read bond list from dump file
d.extra(data) extract bond/tri/line list from data
iterator() loops over selected timesteps
iterator() called with arg = 0 first time, with arg = 1 on subsequent calls

View File

@ -14,12 +14,12 @@ from __future__ import print_function
oneline = "Create plots via GnuPlot plotting program"
docstr = """
g = gnu() start up GnuPlot
g.stop() shut down GnuPlot process
g = gnu() start up GnuPlot
g.stop() shut down GnuPlot process
g.plot(a) plot vector A against linear index
g.plot(a,b) plot B against A
g.plot(a,b,c,d,...) plot B against A, D against C, etc
g.plot(a,b) plot B against A
g.plot(a,b,c,d,...) plot B against A, D against C, etc
g.mplot(M,N,S,"file",a,b,...) multiple plots saved to file0000.eps, etc
each plot argument can be a tuple, list, or Numeric/NumPy vector
@ -32,21 +32,21 @@ g.mplot(M,N,S,"file",a,b,...) multiple plots saved to file0000.eps, etc
g("plot 'file.dat' using 2:3 with lines") execute string in GnuPlot
g.enter() enter GnuPlot shell
g.enter() enter GnuPlot shell
gnuplot> plot sin(x) with lines type commands directly to GnuPlot
gnuplot> exit, quit exit GnuPlot shell
gnuplot> exit, quit exit GnuPlot shell
g.export("data",range(100),a,...) create file with columns of numbers
all vectors must be of equal length
could plot from file with GnuPlot command: plot 'data' using 1:2 with lines
g.select(N) figure N becomes the current plot
g.select(N) figure N becomes the current plot
subsequent commands apply to this plot
g.hide(N) delete window for figure N
g.save("file") save current plot as file.eps
g.hide(N) delete window for figure N
g.save("file") save current plot as file.eps
Set attributes for current plot:

View File

@ -6,6 +6,8 @@
# certain rights in this software. This software is distributed under
# the GNU General Public License.
from __future__ import print_function
# log tool
oneline = "Read LAMMPS log files and extract thermodynamic data"
@ -64,7 +66,7 @@ class log:
# --------------------------------------------------------------------
def __init__(self,*list):
def __init__(self,*arglist):
self.nvec = 0
self.names = []
self.ptr = {}
@ -72,18 +74,18 @@ class log:
# flist = list of all log file names
words = list[0].split()
words = arglist[0].split()
self.flist = []
for word in words: self.flist += glob.glob(word)
if len(self.flist) == 0 and len(list) == 1:
raise StandardError,"no log file specified"
if len(self.flist) == 0 and len(arglist) == 1:
raise Exception("no log file specified")
if len(list) == 1:
if len(arglist) == 1:
self.increment = 0
self.read_all()
else:
if len(self.flist) > 1:
raise StandardError,"can only incrementally read one log file"
raise Exception("can only incrementally read one log file")
self.increment = 1
self.eof = 0
@ -92,24 +94,23 @@ class log:
def read_all(self):
self.read_header(self.flist[0])
if self.nvec == 0: raise StandardError,"log file has no values"
if self.nvec == 0: raise Exception("log file has no values")
# read all files
for file in self.flist: self.read_one(file)
print
# sort entries by timestep, cull duplicates
self.data.sort(self.compare)
self.data.sort(key=(lambda elem: elem[0]))
self.cull()
self.nlen = len(self.data)
print "read %d log entries" % self.nlen
print("read %d log entries" % self.nlen)
# --------------------------------------------------------------------
def next(self):
if not self.increment: raise StandardError,"cannot read incrementally"
if not self.increment: raise Exception("cannot read incrementally")
if self.nvec == 0:
try: open(self.flist[0],'r')
@ -124,12 +125,12 @@ class log:
def get(self,*keys):
if len(keys) == 0:
raise StandardError, "no log vectors specified"
raise Exception( "no log vectors specified" )
map = []
colmap = []
for key in keys:
if self.ptr.has_key(key):
map.append(self.ptr[key])
if key in self.ptr:
colmap.append(self.ptr[key])
else:
count = 0
for i in range(self.nvec):
@ -137,27 +138,27 @@ class log:
count += 1
index = i
if count == 1:
map.append(index)
colmap.append(index)
else:
raise StandardError, "unique log vector %s not found" % key
raise Exception( "unique log vector %s not found" % key)
vecs = []
for i in range(len(keys)):
vecs.append(self.nlen * [0])
for j in xrange(self.nlen):
vecs[i][j] = self.data[j][map[i]]
for j in range(self.nlen):
vecs[i][j] = self.data[j][colmap[i]]
if len(keys) == 1: return vecs[0]
else: return vecs
# --------------------------------------------------------------------
def write(self,filename,*keys):
def write(self,filename,writenames,*keys):
if len(keys):
map = []
colmap = []
for key in keys:
if self.ptr.has_key(key):
map.append(self.ptr[key])
if key in self.ptr:
colmap.append(self.ptr[key])
else:
count = 0
for i in range(self.nvec):
@ -165,17 +166,27 @@ class log:
count += 1
index = i
if count == 1:
map.append(index)
colmap.append(index)
else:
raise StandardError, "unique log vector %s not found" % key
raise Exception( "unique log vector %s not found" % key)
else:
map = range(self.nvec)
colmap = range(self.nvec)
f = open(filename,"w")
for i in xrange(self.nlen):
for j in xrange(len(map)):
print >>f,self.data[i][map[j]],
print >>f
# write col names from dict in the right order
if writenames:
print("# ", file=f, end="")
colnames = [k for j in colmap for k,v in self.ptr.items() if v == j]
for j in range(len(colnames)):
print(colnames[j], file=f, end=" ")
print("\n", file=f, end="")
# write data
for i in range(self.nlen):
for j in range(len(colmap)):
print(self.data[i][colmap[j]],file=f,end=" "),
print("\n",file=f,end="")
f.close()
# --------------------------------------------------------------------
@ -243,18 +254,18 @@ class log:
# --------------------------------------------------------------------
def read_one(self,*list):
def read_one(self,*arglist):
# if 2nd arg exists set file ptr to that value
# if 2nd arg exists set file ptr io that value
# read entire (rest of) file into txt
file = list[0]
file = arglist[0]
if file[-3:] == ".gz":
f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'rb')
else:
f = open(file,'rb')
if len(list) == 2: f.seek(list[1])
if len(arglist) == 2: f.seek(arglist[1])
txt = f.read()
if file[-3:] == ".gz": eof = 0
else: eof = f.tell()
@ -270,12 +281,12 @@ class log:
# set start = position in file to start looking for next chunk
# rewind eof if final entry is incomplete
s1 = txt.find(self.firststr,start)
s2 = txt.find("Loop time of",start+1)
s1 = txt.find(self.firststr.encode('utf-8'),start)
s2 = txt.find("Loop time of".encode('utf-8'),start+1)
if s1 >= 0 and s2 >= 0 and s1 < s2: # found s1,s2 with s1 before s2
if self.style == 2:
s1 = txt.find("\n",s1) + 1
s1 = txt.find("\n".encode('utf-8'),s1) + 1
elif s1 >= 0 and s2 >= 0 and s2 < s1: # found s1,s2 with s2 before s1
s1 = 0
elif s1 == -1 and s2 >= 0: # found s2, but no s1
@ -284,25 +295,25 @@ class log:
elif s1 >= 0 and s2 == -1: # found s1, but no s2
last = 1
if self.style == 1:
s2 = txt.rfind("\n--",s1) + 1
s2 = txt.rfind("\n--".encode('utf-8'),s1) + 1
else:
s1 = txt.find("\n",s1) + 1
s2 = txt.rfind("\n",s1) + 1
s1 = txt.find("\n".encode('utf-8'),s1) + 1
s2 = txt.rfind("\n".encode('utf-8'),s1) + 1
eof -= len(txt) - s2
elif s1 == -1 and s2 == -1: # found neither
# could be end-of-file section
# or entire read was one chunk
if txt.find("Loop time of",start) == start: # end of file, so exit
if txt.find("Loop time of".encode('utf-8'),start) == start: # end of file, so exit
eof -= len(txt) - start # reset eof to "Loop"
break
last = 1 # entire read is a chunk
s1 = 0
if self.style == 1:
s2 = txt.rfind("\n--",s1) + 1
s2 = txt.rfind("\n--".encode('utf-8'),s1) + 1
else:
s2 = txt.rfind("\n",s1) + 1
s2 = txt.rfind("\n".encode('utf-8'),s1) + 1
eof -= len(txt) - s2
if s1 == s2: break
@ -313,24 +324,23 @@ class log:
# parse each entry for numeric fields, append to data
if self.style == 1:
sections = chunk.split("\n--")
sections = chunk.split("\n--".encode('utf-8'))
pat1 = re.compile("Step\s*(\S*)\s")
pat2 = re.compile("=\s*(\S*)")
for section in sections:
word1 = [re.search(pat1,section).group(1)]
word2 = re.findall(pat2,section)
words = word1 + word2
self.data.append(map(float,words))
self.data.append(list(map(float,words)))
else:
lines = chunk.split("\n")
lines = chunk.split("\n".encode('utf-8'))
for line in lines:
words = line.split()
self.data.append(map(float,words))
self.data.append(list(map(float,words)))
# print last timestep of chunk
print int(self.data[len(self.data)-1][0]),
print(int(self.data[len(self.data)-1][0]),)
sys.stdout.flush()
return eof

View File

@ -6,6 +6,8 @@
# certain rights in this software. This software is distributed under
# the GNU General Public License.
from __future__ import print_function
# xyz tool
oneline = "Convert LAMMPS snapshots to XYZ format"
@ -56,18 +58,18 @@ class xyz:
if flag == -1: break
time,box,atoms,bonds,tris,lines = self.data.viz(which)
print >>f,len(atoms)
print >>f,"Atoms"
print(len(atoms), file=f)
print("Atoms", file=f)
for atom in atoms:
itype = int(atom[1])
print >>f,itype,atom[2],atom[3],atom[4]
print(itype,atom[2],atom[3],atom[4], file=f)
print time,
print(time)
sys.stdout.flush()
n += 1
f.close()
print "\nwrote %d snapshots to %s in XYZ format" % (n,file)
print("\nwrote %d snapshots to %s in XYZ format" % (n,file))
# --------------------------------------------------------------------
@ -91,17 +93,17 @@ class xyz:
file = root + str(n)
file += ".xyz"
f = open(file,"w")
print >>f,len(atoms)
print >>f,"Atoms"
print(len(atoms), file=f)
print("Atoms", file=f)
for atom in atoms:
itype = int(atom[1])
print >>f,itype,atom[2],atom[3],atom[4]
print time,
print(itype,atom[2],atom[3],atom[4], file=f)
print(time)
sys.stdout.flush()
f.close()
n += 1
print "\nwrote %s snapshots in XYZ format" % n
print("\nwrote %s snapshots in XYZ format" % n)
# --------------------------------------------------------------------
@ -113,9 +115,9 @@ class xyz:
which = self.data.findtime(time)
time,box,atoms,bonds,tris,lines = self.data.viz(which)
f = open(file,"w")
print >>f,len(atoms)
print >>f,"Atoms"
print(len(atoms), file=f)
print("Atoms", file=f)
for atom in atoms:
itype = int(atom[1])
print >>f,itype,atom[2],atom[3],atom[4]
print(itype,atom[2],atom[3],atom[4], file=f)
f.close()