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

Merge pull request #4183 from alphataubio/spica-kk

Browse Source 
Add Kokkos version of spica pair and angle styles
This commit is contained in:
Axel Kohlmeyer
2024-07-26 22:45:47 -04:00
committed by GitHub
parent 1014cd1e0e 2a99b9957d
commit 7396426963
29 changed files with 2637 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_bond.rst
View File

@ -101,7 +101,7 @@ OPT.
* :doc:`mesocnt <angle_mesocnt>` * :doc:`mesocnt <angle_mesocnt>`
* :doc:`mm3 <angle_mm3>` * :doc:`mm3 <angle_mm3>`
* :doc:`quartic (o) <angle_quartic>` * :doc:`quartic (o) <angle_quartic>`
* :doc:`spica (o) <angle_spica>` * :doc:`spica (ko) <angle_spica>`
* :doc:`table (o) <angle_table>` * :doc:`table (o) <angle_table>`
.. _dihedral: .. _dihedral:

2
doc/src/Commands_pair.rst
View File

@ -195,7 +195,7 @@ OPT.
* :doc:`lj/mdf <pair_mdf>` * :doc:`lj/mdf <pair_mdf>`
* :doc:`lj/relres (o) <pair_lj_relres>` * :doc:`lj/relres (o) <pair_lj_relres>`
* :doc:`lj/spica (gko) <pair_spica>` * :doc:`lj/spica (gko) <pair_spica>`
* :doc:`lj/spica/coul/long (go) <pair_spica>` * :doc:`lj/spica/coul/long (gko) <pair_spica>`
* :doc:`lj/spica/coul/msm (o) <pair_spica>` * :doc:`lj/spica/coul/msm (o) <pair_spica>`
* :doc:`lj/sf/dipole/sf (go) <pair_dipole>` * :doc:`lj/sf/dipole/sf (go) <pair_dipole>`
* :doc:`lj/smooth (go) <pair_lj_smooth>` * :doc:`lj/smooth (go) <pair_lj_smooth>`

3
doc/src/angle_spica.rst
View File

@ -1,10 +1,11 @@
.. index:: angle_style spica .. index:: angle_style spica
.. index:: angle_style spica/omp .. index:: angle_style spica/omp
.. index:: angle_style spica/kk
angle_style spica command angle_style spica command
========================= =========================
Accelerator Variants: *spica/omp* Accelerator Variants: *spica/omp*, *spica/kk*
Syntax Syntax
"""""" """"""

3
doc/src/pair_spica.rst
View File

@ -5,6 +5,7 @@
.. index:: pair_style lj/spica/coul/long .. index:: pair_style lj/spica/coul/long
.. index:: pair_style lj/spica/coul/long/gpu .. index:: pair_style lj/spica/coul/long/gpu
.. index:: pair_style lj/spica/coul/long/omp .. index:: pair_style lj/spica/coul/long/omp
.. index:: pair_style lj/spica/coul/long/kk
.. index:: pair_style lj/spica/coul/msm .. index:: pair_style lj/spica/coul/msm
.. index:: pair_style lj/spica/coul/msm/omp .. index:: pair_style lj/spica/coul/msm/omp
@ -16,7 +17,7 @@ Accelerator Variants: *lj/spica/gpu*, *lj/spica/kk*, *lj/spica/omp*
pair_style lj/spica/coul/long command pair_style lj/spica/coul/long command
===================================== =====================================
Accelerator Variants: *lj/spica/coul/long/gpu*, *lj/spica/coul/long/omp* Accelerator Variants: *lj/spica/coul/long/gpu*, *lj/spica/coul/long/omp*, *lj/spica/coul/long/kk*
pair_style lj/spica/coul/msm command pair_style lj/spica/coul/msm command
==================================== ====================================

4
examples/PACKAGES/cgspica/sds-monolayer/in.sds-regular
View File

@ -5,9 +5,9 @@ dimension 3
atom_style full atom_style full
processors * * 1 processors * * 1
pair_style lj/sdk/coul/long 15.0 # compatible with "lj/spica/coul/long" pair_style lj/spica/coul/long 15.0
bond_style harmonic bond_style harmonic
angle_style sdk # compatible with "spica" angle_style spica
special_bonds lj/coul 0.0 0.0 1.0 special_bonds lj/coul 0.0 0.0 1.0
read_data data.sds.gz read_data data.sds.gz

340
examples/charmmfsw/charmmff.cmap Normal file
View File

@ -0,0 +1,340 @@
#CMAP for C NH1 CT1 C NH1 CT1 C NH1; id=1
#phi = -180.000000
0.130000 0.770000 0.970000 1.250000 2.120000
2.720000 2.090000 1.790000 0.780000 -0.690000
1.000000 -2.200000 -4.830000 -4.820000 -4.910000
-3.590000 -2.770000 -2.780000 -2.450000 -2.350000
-2.340000 -1.520000 -0.950000 -0.040000
#phi = -165.000000
-0.130000 1.380000 1.580000 1.870000 2.400000
2.490000 2.440000 1.930000 1.090000 0.640000
0.260000 -2.800000 -4.010000 -4.140000 -3.420000
-2.600000 -2.300000 -1.500000 -1.100000 -0.860000
-0.640000 -0.210000 -1.080000 -1.120000
#phi = -150.000000
0.080000 1.420000 1.620000 2.050000 2.650000
2.720000 2.320000 1.990000 1.560000 2.460000
-0.230000 -1.820000 -2.580000 -3.010000 -2.550000
-1.890000 -1.350000 -0.730000 0.070000 -0.230000
-0.770000 -1.280000 -1.290000 -0.820000
#phi = -135.000000
0.930000 1.520000 2.240000 2.550000 3.110000
2.920000 2.460000 2.190000 2.060000 1.850000
0.120000 -1.180000 -2.000000 -2.280000 -1.960000
-1.340000 -0.930000 0.020000 0.310000 -0.520000
-1.150000 -0.980000 -0.570000 -0.440000
#phi = -120.000000
1.360000 1.960000 2.700000 3.040000 3.700000
3.560000 2.640000 2.770000 2.720000 1.630000
0.710000 -0.790000 -2.120000 -2.630000 -1.800000
-0.430000 -0.060000 0.440000 0.910000 -0.550000
-0.970000 -0.860000 -0.250000 0.450000
#phi = -105.000000
2.050000 2.540000 2.820000 3.090000 3.370000
3.550000 3.070000 2.900000 2.960000 2.120000
0.910000 -0.820000 -2.090000 -2.240000 -1.460000
0.210000 0.080000 0.770000 1.040000 -0.120000
-0.320000 -0.160000 0.310000 0.730000
#phi = -90.000000
1.450000 2.750000 2.740000 3.160000 3.450000
3.340000 3.180000 3.900000 3.340000 2.440000
0.910000 -0.610000 -1.510000 -1.620000 -0.960000
-0.020000 0.420000 0.910000 0.460000 0.150000
-0.070000 0.020000 0.280000 0.750000
#phi = -75.000000
1.380000 3.350000 2.350000 3.060000 3.810000
3.700000 3.580000 4.210000 3.540000 1.690000
0.100000 -0.680000 -0.120000 -0.430000 -0.600000
0.230000 0.420000 0.300000 0.550000 0.190000
-0.250000 -0.190000 -0.250000 0.470000
#phi = -60.000000
0.240000 1.230000 1.720000 3.170000 4.210000
4.390000 4.280000 3.670000 2.270000 -0.480000
-0.410000 -0.040000 -0.360000 -0.820000 -0.170000
0.140000 0.270000 0.320000 0.310000 -0.670000
-0.950000 -1.530000 -1.480000 -0.200000
#phi = -45.000000
-1.180000 0.080000 2.350000 4.210000 5.380000
5.390000 4.380000 2.460000 1.120000 0.110000
0.010000 -0.150000 -0.800000 -0.580000 0.080000
0.270000 -0.050000 0.380000 0.250000 -0.890000
-1.580000 -1.950000 -1.980000 -2.000000
#phi = -30.000000
-1.170000 1.070000 4.180000 6.740000 6.070000
4.810000 2.780000 1.320000 0.770000 -0.010000
0.280000 -0.710000 1.310000 1.520000 1.920000
2.220000 0.190000 0.530000 0.330000 -1.600000
-2.850000 -3.550000 -3.280000 -2.660000
#phi = -15.000000
0.290000 5.590000 3.730000 3.220000 3.270000
2.520000 1.590000 1.380000 0.860000 0.660000
1.620000 0.850000 0.510000 0.740000 1.020000
1.620000 -0.340000 0.180000 -0.610000 -2.560000
-3.790000 -3.810000 -3.160000 -1.750000
#phi = 0.000000
2.830000 0.790000 0.320000 0.480000 0.630000
0.980000 1.240000 1.670000 1.650000 2.520000
1.610000 0.780000 0.120000 0.070000 0.120000
-1.570000 -1.210000 -1.930000 -2.600000 -3.790000
-3.930000 -3.620000 -2.680000 -0.920000
#phi = 15.000000
-0.780000 -1.910000 -2.050000 -1.850000 -1.050000
0.180000 1.680000 2.220000 1.360000 2.450000
1.440000 0.680000 -0.240000 -0.540000 -0.790000
-2.180000 -3.210000 -4.350000 -3.940000 -3.910000
-3.460000 -2.770000 1.760000 0.310000
#phi = 30.000000
-2.960000 -3.480000 -3.440000 -2.400000 -1.130000
0.340000 1.430000 1.390000 0.970000 2.460000
1.520000 0.550000 -0.410000 -1.480000 -3.580000
-4.130000 -4.560000 -4.440000 -3.580000 -2.960000
-1.960000 -1.070000 -1.600000 -2.450000
#phi = 45.000000
-4.020000 -3.840000 -3.370000 -2.330000 -0.980000
0.360000 0.810000 0.750000 0.500000 1.900000
0.770000 -0.420000 -3.290000 -3.910000 -4.520000
-4.890000 -3.850000 -4.150000 -2.670000 -2.370000
-2.860000 -3.420000 -3.670000 -3.600000
#phi = 60.000000
-3.350000 -2.980000 -2.320000 -1.240000 -0.260000
0.720000 0.670000 0.440000 2.400000 1.630000
-2.010000 -3.310000 -3.990000 -4.530000 -4.850000
-3.770000 -3.940000 -3.890000 -2.610000 -3.510000
-3.760000 -3.640000 -3.450000 -3.340000
#phi = 75.000000
-2.250000 -1.640000 -1.010000 0.040000 0.640000
0.820000 0.520000 -0.010000 -0.370000 -1.190000
-2.390000 -3.380000 -4.500000 -5.590000 -5.510000
-4.940000 -3.830000 -3.840000 -3.700000 -4.150000
-4.170000 -3.730000 -3.740000 -2.620000
#phi = 90.000000
-1.720000 -1.180000 -0.430000 0.280000 0.810000
0.800000 0.480000 -0.340000 -0.790000 -1.770000
-2.810000 -3.800000 -5.220000 -6.280000 -6.580000
-5.640000 -5.060000 -4.020000 -4.150000 -4.470000
-4.100000 -3.770000 -3.160000 -2.650000
#phi = 105.000000
-1.850000 -1.090000 -0.450000 0.130000 1.010000
0.880000 0.490000 -0.220000 -0.860000 -1.680000
-3.010000 -4.130000 -5.990000 -6.860000 -6.830000
-5.850000 -3.860000 -4.860000 -4.910000 -4.720000
-4.600000 -4.090000 -3.270000 -2.410000
#phi = 120.000000
-1.970000 -1.120000 -0.540000 -0.150000 0.760000
1.040000 0.760000 0.310000 -0.330000 -1.870000
-3.370000 -5.010000 -6.120000 -7.050000 -6.980000
-3.700000 -4.510000 -5.090000 -5.420000 -4.850000
-4.440000 -4.000000 -3.420000 -2.750000
#phi = 135.000000
-2.110000 -1.170000 -0.320000 -0.010000 0.320000
1.090000 0.940000 0.630000 -0.170000 -1.830000
-3.470000 -4.950000 -6.110000 -1.920000 -4.050000
-5.000000 -5.000000 -4.840000 -4.890000 -4.300000
-4.490000 -4.440000 -4.160000 -3.180000
#phi = 150.000000
-1.760000 -0.400000 0.020000 0.360000 0.630000
1.260000 1.360000 0.950000 -0.070000 -1.480000
-3.150000 1.840000 -1.760000 -5.090000 -5.740000
-5.390000 -4.780000 -4.190000 -4.120000 -4.040000
-4.130000 -4.030000 -4.030000 -2.940000
#phi = 165.000000
-0.810000 -0.070000 0.380000 0.540000 1.280000
1.640000 1.700000 1.520000 0.630000 -1.090000
-2.740000 -0.740000 -4.560000 -6.410000 -5.890000
-5.140000 -4.190000 -3.670000 -3.840000 -3.560000
-3.550000 -3.250000 -2.750000 -1.810000
#CMAP for C NH1 CT2 C NH1 CT2 C NH1; id=2
#phi = -180.000000
0.235350 0.182300 0.177200 0.396800 0.859400
1.489700 2.092500 2.297700 1.808600 0.696200
-0.563300 -1.432700 -1.015100 1.426300 -0.564300
0.696200 1.808200 2.301700 2.092600 1.489100
0.859500 0.396900 0.176900 0.182400
#phi = -165.000000
0.020100 -0.203800 -0.269700 0.014200 0.620800
1.392400 2.046200 2.188200 1.683900 0.688500
-0.373700 -0.703500 0.837800 3.704000 -0.730100
0.594100 1.713100 2.205800 2.026400 1.529800
1.027400 0.623800 0.348400 0.182800
#phi = -150.000000
-0.533600 -0.807400 -0.804600 -0.379800 0.365300
1.168000 1.641000 1.618100 1.302200 0.615100
0.065700 0.738500 2.959500 -2.036600 -0.934600
0.407900 1.517000 1.984800 1.833100 1.435200
0.995600 0.562200 0.150600 -0.209000
#phi = -135.000000
-1.208500 -1.429400 -1.319200 -0.817500 -0.112400
0.454400 0.737600 0.879300 0.850100 0.670300
0.943500 -2.651200 -2.829400 -2.199100 -1.065700
0.279600 1.322000 1.668300 1.521300 1.193900
0.765300 0.246000 -0.315500 -0.823200
#phi = -120.000000
-1.789100 -1.965500 -1.860700 -1.447900 -0.896500
-0.401000 -0.015100 0.321300 0.634600 0.976300
-1.977500 -2.883200 -2.848500 -2.137900 -0.960300
0.308700 1.098100 1.245300 1.133600 0.881800
0.448200 -0.153900 -0.823700 -1.404300
#phi = -105.000000
-2.246700 -2.487000 -2.473700 -2.135600 -1.577700
-0.980600 -0.429100 0.144700 0.734000 -0.918300
-2.299200 -2.882200 -2.668600 -1.847100 -0.719800
0.107000 0.496000 0.553500 0.584300 0.494000
0.098300 -0.529800 -1.237900 -1.840100
#phi = -90.000000
-2.851100 -3.181100 -3.199500 -2.785300 -2.054300
-1.242900 -0.476500 0.288100 -0.045300 -1.470600
-2.558800 -2.869400 -2.450300 -1.582200 -0.930800
-0.426400 -0.022700 0.000000 -0.097400 -0.136100
-0.439600 -1.038600 -1.741000 -2.373200
#phi = -75.000000
-3.961800 -4.268200 -4.109000 -3.364700 -2.252200
-1.140400 -0.209800 0.487300 -0.746200 -2.127700
-2.932100 -2.898500 -2.247900 -1.730400 -1.177200
-0.448200 0.034900 -0.073300 -0.531600 -0.933300
-1.360700 -2.009200 -2.745700 -3.424900
#phi = -60.000000
-5.408000 -5.355100 -4.640100 -3.283200 -1.710200
-0.423800 0.354400 -0.103700 -1.577700 -2.828300
-3.151200 -2.649200 -2.183000 -1.761200 -0.981700
-0.174700 0.262600 0.039200 -0.663000 -1.530700
-2.478200 -3.465600 -4.334200 -5.011200
#phi = -45.000000
-6.093200 -5.298400 -3.816620 -1.922530 -0.196160
0.768200 0.568500 -0.831300 -2.343900 -3.037100
-2.663700 -2.191100 -2.022900 -1.438500 -0.649000
0.077000 0.441500 0.257500 -0.491100 -1.820600
-3.473100 -4.895200 -5.790700 -6.205900
#phi = -30.000000
-5.258225 -3.675795 -1.631110 0.430085 1.496470
0.318200 -0.555100 -1.695500 -2.434200 -2.192600
-1.691300 -1.890000 -1.708500 -1.206300 -0.567400
0.054300 0.497200 0.599600 -0.171000 -2.137600
-4.237000 -5.584100 -6.135100 -6.067000
#phi = -15.000000
-3.161820 -0.902080 1.432450 -1.452885 -1.560780
-1.665600 -1.783100 -1.755100 -1.329300 -0.731100
-1.317000 -1.662800 -1.601200 -1.294900 -0.817300
-0.197100 0.549500 0.850400 -0.689700 -2.819900
-4.393000 -5.111500 -5.205690 -4.654785
#phi = 0.000000
0.034035 -2.349860 -3.412065 -3.620070 -3.450950
-2.875650 -1.787800 -0.541250 0.410450 -0.372500
-1.126850 -1.498450 -1.608700 -1.498450 -1.126850
-0.372500 0.410450 -0.541250 -1.787800 -2.875650
-3.450950 -3.620070 -3.412065 -2.349860
#phi = 15.000000
-3.162345 -4.654785 -5.205690 -5.111500 -4.393000
-2.819900 -0.689700 0.850400 0.549500 -0.197100
-0.817300 -1.294900 -1.601200 -1.662800 -1.317000
-0.731100 -1.329300 -1.755100 -1.783100 -1.665600
-1.560780 -1.452885 1.432450 -0.902080
#phi = 30.000000
-5.258220 -6.067000 -6.135100 -5.584100 -4.237000
-2.137600 -0.171000 0.599600 0.497200 0.054300
-0.567400 -1.206300 -1.708500 -1.890000 -1.691300
-2.192600 -2.434200 -1.695500 -0.555100 0.318200
1.496470 0.430085 -1.631110 -3.675795
#phi = 45.000000
-6.093300 -6.205900 -5.790700 -4.895200 -3.473100
-1.820600 -0.491100 0.257500 0.441500 0.077000
-0.649000 -1.438500 -2.022900 -2.191100 -2.663700
-3.037100 -2.343900 -0.831300 0.568500 0.768200
-0.196160 -1.922530 -3.816620 -5.298400
#phi = 60.000000
-5.407500 -5.011200 -4.334200 -3.465600 -2.478200
-1.530700 -0.663000 0.039200 0.262600 -0.174700
-0.981700 -1.761200 -2.183000 -2.649200 -3.151200
-2.828300 -1.577700 -0.103700 0.354400 -0.423800
-1.710200 -3.283200 -4.640100 -5.355100
#phi = 75.000000
-3.961900 -3.424900 -2.745700 -2.009200 -1.360700
-0.933300 -0.531600 -0.073300 0.034900 -0.448200
-1.177200 -1.730400 -2.247900 -2.898500 -2.932100
-2.127700 -0.746200 0.487300 -0.209800 -1.140400
-2.252200 -3.364700 -4.109000 -4.268200
#phi = 90.000000
-2.854500 -2.373200 -1.741000 -1.038600 -0.439600
-0.136100 -0.097400 0.000000 -0.022700 -0.426400
-0.930800 -1.582200 -2.450300 -2.869400 -2.558800
-1.470600 -0.045300 0.288100 -0.476500 -1.242900
-2.054300 -2.785300 -3.199500 -3.181100
#phi = 105.000000
-2.246400 -1.840100 -1.237900 -0.529800 0.098300
0.494000 0.584300 0.553500 0.496000 0.107000
-0.719800 -1.847100 -2.668600 -2.882200 -2.299200
-0.918300 0.734000 0.144700 -0.429100 -0.980600
-1.577700 -2.135600 -2.473700 -2.487000
#phi = 120.000000
-1.788800 -1.404300 -0.823700 -0.153900 0.448200
0.881800 1.133600 1.245300 1.098100 0.308700
-0.960300 -2.137900 -2.848500 -2.883200 -1.977500
0.976300 0.634600 0.321300 -0.015100 -0.401000
-0.896500 -1.447900 -1.860700 -1.965500
#phi = 135.000000
-1.208900 -0.823200 -0.315500 0.246000 0.765300
1.193900 1.521300 1.668300 1.322000 0.279600
-1.065700 -2.199100 -2.829400 -2.651200 0.943500
0.670300 0.850100 0.879300 0.737600 0.454400
-0.112400 -0.817500 -1.319200 -1.429400
#phi = 150.000000
-0.533400 -0.209000 0.150600 0.562200 0.995600
1.435200 1.833100 1.984800 1.517000 0.407900
-0.934600 -2.036600 2.959500 0.738500 0.065700
0.615100 1.302200 1.618100 1.641000 1.168000
0.365300 -0.379800 -0.804600 -0.807400
#phi = 165.000000
0.019900 0.182800 0.348400 0.623800 1.027400
1.529800 2.026400 2.205800 1.713100 0.594100
-0.730100 3.704000 0.837800 -0.703500 -0.373700
0.688500 1.683900 2.188200 2.046200 1.392400
0.620800 0.014200 -0.269700 -0.203800

BIN
examples/charmmfsw/data.charmmfsw.gz Normal file
View File

Binary file not shown.

37
examples/charmmfsw/in.charmmfsw Normal file
View File

@ -0,0 +1,37 @@
# charmmfsw example generated by https://charmm-gui.org/
# from PDB structure 1HVN (https://www.rcsb.org/structure/1HVN)
#
# Dependencies: packages MOLECULE / KSPACE / RIGID
# To test with KOKKOS: lmp -k on g 1 -sf kk -pk kokkos neigh half -in in.charmmfsw
units real
boundary p p p
newton off
pair_style lj/charmmfsw/coul/long 10 12
pair_modify mix arithmetic
kspace_style pppm 1e-6
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
special_bonds charmm
improper_style harmonic
timestep 2
fix cmap all cmap charmmff.cmap
fix_modify cmap energy yes
read_data data.charmmfsw.gz fix cmap crossterm CMAP
neighbor 2 bin
neigh_modify delay 2 every 1
fix 1 all shake 1e-6 100 100 m 1.008 a 142
fix 2 all nvt temp 303.15 303.15 100.0
thermo 10
thermo_style custom step etotal evdwl ecoul elong edihed pe ke temp press
run 100

221
examples/charmmfsw/log.26Jul24.charmmfsw.g++.1 Normal file
View File

@ -0,0 +1,221 @@
LAMMPS (27 Jun 2024)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# charmmfsw example generated by https://charmm-gui.org/
# from PDB structure 1HVN (https://www.rcsb.org/structure/1HVN)
#
# Dependencies: packages MOLECULE / KSPACE / RIGID
# To test with KOKKOS: lmp -k on g 1 -sf kk -pk kokkos neigh half -in in.charmmfsw
units real
boundary p p p
newton off
pair_style lj/charmmfsw/coul/long 10 12
Switching to CHARMM coulomb energy conversion constant (src/KSPACE/pair_lj_charmmfsw_coul_long.cpp:63)
pair_modify mix arithmetic
kspace_style pppm 1e-6
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
special_bonds charmm
improper_style harmonic
timestep 2
fix cmap all cmap charmmff.cmap
Reading CMAP parameters from: charmmff.cmap
Read in CMAP data for 2 crossterm types
fix_modify cmap energy yes
read_data data.charmmfsw.gz fix cmap crossterm CMAP
Reading data file ...
orthogonal box = (-24 -24 -24) to (24 24 24)
1 by 1 by 1 MPI processor grid
reading atoms ...
10245 atoms
reading velocities ...
10245 velocities
scanning bonds ...
4 = max bonds/atom
scanning angles ...
15 = max angles/atom
scanning dihedrals ...
48 = max dihedrals/atom
scanning impropers ...
4 = max impropers/atom
orthogonal box = (-24 -24 -24) to (24 24 24)
1 by 1 by 1 MPI processor grid
reading bonds ...
6973 bonds
reading angles ...
4057 angles
reading dihedrals ...
1363 dihedrals
reading impropers ...
70 impropers
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
4 = max # of 1-2 neighbors
9 = max # of 1-3 neighbors
17 = max # of 1-4 neighbors
21 = max # of special neighbors
special bonds CPU = 0.002 seconds
read_data CPU = 0.072 seconds
neighbor 2 bin
neigh_modify delay 2 every 1
fix 1 all shake 1e-6 100 100 m 1.008 a 142
Finding SHAKE clusters ...
75 = # of size 2 clusters
47 = # of size 3 clusters
9 = # of size 4 clusters
3265 = # of frozen angles
find clusters CPU = 0.003 seconds
fix 2 all nvt temp 303.15 303.15 100.0
thermo 10
thermo_style custom step etotal evdwl ecoul elong edihed pe ke temp press
run 100
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
@Article{Gissinger24,
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
journal = {J. Phys. Chem. B},
year = 2024,
volume = 128,
number = 13,
pages = {3282–-3297}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.27938162
grid = 54 54 54
stencil order = 5
estimated absolute RMS force accuracy = 0.00036407395
estimated relative force accuracy = 1.0963718e-06
using double precision FFTW3
3d grid and FFT values/proc = 226981 157464
Generated 2016 of 2016 mixed pair_coeff terms from arithmetic mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 2 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/charmmfsw/coul/long, perpetual
attributes: half, newton off
pair build: half/bin/newtoff
stencil: full/bin/3d
bin: standard
SHAKE stats (type/ave/delta/count) on step 0
Bond: 16 1.09 1.38032e-07 6
Bond: 18 1.09 1.00046e-07 3
Bond: 34 1.111 1.11388e-06 10
Bond: 39 1.111 4.83041e-08 5
Bond: 43 1.111 1.97842e-07 10
Bond: 44 1.111 1.71815e-06 10
Bond: 59 1.111 8.42509e-08 2
Bond: 62 1.111 2.84854e-08 2
Bond: 63 1.111 2.14153e-07 46
Bond: 64 1.111 1.59305e-07 18
Bond: 65 1.08 5.67061e-07 16
Bond: 66 1.08 1.43965e-06 4
Bond: 67 1 1.81926e-07 10
Bond: 68 1.01 0 1
Bond: 69 1.08 1.34571e-07 5
Bond: 70 1.09 0 1
Bond: 71 1.083 0 1
Bond: 72 0.9572 2.71955e-07 6530
Bond: 75 1 1.46045e-07 10
Bond: 79 0.997 5.24499e-07 17
Bond: 81 1 1.32984e-07 4
Bond: 84 1.04 7.65389e-07 9
Bond: 87 1 0 1
Bond: 95 0.96 5.75241e-07 2
Bond: 97 1.325 4.3613e-08 3
Angle: 142 104.52 2.67611e-05 3265
Per MPI rank memory allocation (min/avg/max) = 143.6 | 143.6 | 143.6 Mbytes
Step TotEng E_vdwl E_coul E_long E_dihed PotEng KinEng Temp Press
0 -27877.652 3447.5013 144035.68 -182420.51 343.05623 -34213.5 6335.8474 307.44113 -989.27065
10 -27879.086 3334.4154 144205.4 -182416.19 348.14696 -34133.566 6254.4808 303.49289 -1211.2863
20 -27882.193 3293.7931 144272.04 -182415.87 333.20456 -34116.91 6234.7164 302.53384 -1041.5231
30 -27886.779 3177.7183 144344.61 -182409.28 340.77044 -34166.241 6279.462 304.70508 -1538.0247
40 -27892.698 3186.4294 144409.85 -182417.01 337.80177 -34097.62 6204.9214 301.08807 -1516.1201
50 -27898.215 3198.5531 144426.3 -182405.24 336.58074 -34049.909 6151.6947 298.50529 -1349.4431
60 -27900.589 3163.4592 144400.32 -182414.85 341.17705 -34110.926 6210.3369 301.35085 -1695.3697
70 -27900.487 3223.7183 144242.71 -182409.21 341.09496 -34188.493 6288.0059 305.11967 -1493.2031
80 -27901.07 3274.244 144265.07 -182417.68 344.0409 -34177.343 6276.2725 304.55032 -1273.0263
90 -27905.672 3237.6056 144288.71 -182418.22 342.15013 -34187.814 6282.1417 304.83511 -1268.0436
SHAKE stats (type/ave/delta/count) on step 100
Bond: 16 1.09 3.78281e-07 6
Bond: 18 1.09 1.12288e-07 3
Bond: 34 1.111 7.60709e-07 10
Bond: 39 1.111 2.37855e-07 5
Bond: 43 1.111 6.00872e-07 10
Bond: 44 1.111 3.75324e-07 10
Bond: 59 1.111 1.12311e-07 2
Bond: 62 1.111 2.99471e-07 2
Bond: 63 1.111 6.10589e-07 46
Bond: 64 1.111 4.50733e-07 18
Bond: 65 1.08 2.90668e-07 16
Bond: 66 1.08 1.61592e-07 4
Bond: 67 1 5.4508e-07 10
Bond: 68 1.01 0 1
Bond: 69 1.08 4.1398e-07 5
Bond: 70 1.09 0 1
Bond: 71 1.083 0 1
Bond: 72 0.9572 1.76706e-06 6530
Bond: 75 1 3.96686e-07 10
Bond: 79 0.997 7.72922e-07 17
Bond: 81 1 1.30673e-07 4
Bond: 84 1.04 1.44551e-07 9
Bond: 87 1 0 1
Bond: 95 0.96 1.03526e-07 2
Bond: 97 1.325 3.64689e-08 3
Angle: 142 104.52 0.000130126 3265
100 -27913.241 3159.0677 144299.1 -182414.94 336.48839 -34254.412 6341.1706 307.69943 -1421.2905
Loop time of 11.5304 on 1 procs for 100 steps with 10245 atoms
Performance: 1.499 ns/day, 16.014 hours/ns, 8.673 timesteps/s, 88.852 katom-step/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 | 8.6772 | 8.6772 | 8.6772 | 0.0 | 75.25
Bond | 0.012444 | 0.012444 | 0.012444 | 0.0 | 0.11
Kspace | 1.2286 | 1.2286 | 1.2286 | 0.0 | 10.66
Neigh | 1.5276 | 1.5276 | 1.5276 | 0.0 | 13.25
Comm | 0.010441 | 0.010441 | 0.010441 | 0.0 | 0.09
Output | 0.00055001 | 0.00055001 | 0.00055001 | 0.0 | 0.00
Modify | 0.07101 | 0.07101 | 0.07101 | 0.0 | 0.62
Other | | 0.002628 | | | 0.02
Nlocal: 10245 ave 10245 max 10245 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 30479 ave 30479 max 30479 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 7.05928e+06 ave 7.05928e+06 max 7.05928e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7059275
Ave neighs/atom = 689.04588
Ave special neighs/atom = 2.3664226
Neighbor list builds = 10
Dangerous builds = 0
Total wall time: 0:00:11

221
examples/charmmfsw/log.26Jul24.charmmfsw.g++.4 Normal file
View File

@ -0,0 +1,221 @@
LAMMPS (27 Jun 2024)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# charmmfsw example generated by https://charmm-gui.org/
# from PDB structure 1HVN (https://www.rcsb.org/structure/1HVN)
#
# Dependencies: packages MOLECULE / KSPACE / RIGID
# To test with KOKKOS: lmp -k on g 1 -sf kk -pk kokkos neigh half -in in.charmmfsw
units real
boundary p p p
newton off
pair_style lj/charmmfsw/coul/long 10 12
Switching to CHARMM coulomb energy conversion constant (src/KSPACE/pair_lj_charmmfsw_coul_long.cpp:63)
pair_modify mix arithmetic
kspace_style pppm 1e-6
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
special_bonds charmm
improper_style harmonic
timestep 2
fix cmap all cmap charmmff.cmap
Reading CMAP parameters from: charmmff.cmap
Read in CMAP data for 2 crossterm types
fix_modify cmap energy yes
read_data data.charmmfsw.gz fix cmap crossterm CMAP
Reading data file ...
orthogonal box = (-24 -24 -24) to (24 24 24)
1 by 2 by 2 MPI processor grid
reading atoms ...
10245 atoms
reading velocities ...
10245 velocities
scanning bonds ...
4 = max bonds/atom
scanning angles ...
15 = max angles/atom
scanning dihedrals ...
48 = max dihedrals/atom
scanning impropers ...
4 = max impropers/atom
orthogonal box = (-24 -24 -24) to (24 24 24)
1 by 2 by 2 MPI processor grid
reading bonds ...
6973 bonds
reading angles ...
4057 angles
reading dihedrals ...
1363 dihedrals
reading impropers ...
70 impropers
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
4 = max # of 1-2 neighbors
9 = max # of 1-3 neighbors
17 = max # of 1-4 neighbors
21 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.068 seconds
neighbor 2 bin
neigh_modify delay 2 every 1
fix 1 all shake 1e-6 100 100 m 1.008 a 142
Finding SHAKE clusters ...
75 = # of size 2 clusters
47 = # of size 3 clusters
9 = # of size 4 clusters
3265 = # of frozen angles
find clusters CPU = 0.001 seconds
fix 2 all nvt temp 303.15 303.15 100.0
thermo 10
thermo_style custom step etotal evdwl ecoul elong edihed pe ke temp press
run 100
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
@Article{Gissinger24,
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
journal = {J. Phys. Chem. B},
year = 2024,
volume = 128,
number = 13,
pages = {3282–-3297}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.27938162
grid = 54 54 54
stencil order = 5
estimated absolute RMS force accuracy = 0.00036407395
estimated relative force accuracy = 1.0963718e-06
using double precision FFTW3
3d grid and FFT values/proc = 70516 40824
Generated 2016 of 2016 mixed pair_coeff terms from arithmetic mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 2 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/charmmfsw/coul/long, perpetual
attributes: half, newton off
pair build: half/bin/newtoff
stencil: full/bin/3d
bin: standard
SHAKE stats (type/ave/delta/count) on step 0
Bond: 16 1.09 1.38032e-07 6
Bond: 18 1.09 1.00046e-07 3
Bond: 34 1.111 1.11388e-06 10
Bond: 39 1.111 4.83041e-08 5
Bond: 43 1.111 1.97842e-07 10
Bond: 44 1.111 1.71815e-06 10
Bond: 59 1.111 8.42509e-08 2
Bond: 62 1.111 2.84854e-08 2
Bond: 63 1.111 2.14153e-07 46
Bond: 64 1.111 1.59305e-07 18
Bond: 65 1.08 5.67061e-07 16
Bond: 66 1.08 1.43965e-06 4
Bond: 67 1 1.81926e-07 10
Bond: 68 1.01 0 1
Bond: 69 1.08 1.34571e-07 5
Bond: 70 1.09 0 1
Bond: 71 1.083 0 1
Bond: 72 0.9572 2.71955e-07 6530
Bond: 75 1 1.46045e-07 10
Bond: 79 0.997 5.24499e-07 17
Bond: 81 1 1.32984e-07 4
Bond: 84 1.04 7.65389e-07 9
Bond: 87 1 0 1
Bond: 95 0.96 5.75241e-07 2
Bond: 97 1.325 4.3613e-08 3
Angle: 142 104.52 2.67611e-05 3265
Per MPI rank memory allocation (min/avg/max) = 76.88 | 77.06 | 77.25 Mbytes
Step TotEng E_vdwl E_coul E_long E_dihed PotEng KinEng Temp Press
0 -27877.652 3447.5013 144035.68 -182420.51 343.05623 -34213.5 6335.8474 307.44113 -989.27065
10 -27879.086 3334.4154 144205.4 -182416.19 348.14696 -34133.566 6254.4808 303.49289 -1211.2863
20 -27882.193 3293.7931 144272.04 -182415.87 333.20456 -34116.91 6234.7164 302.53384 -1041.5231
30 -27886.779 3177.7183 144344.61 -182409.28 340.77044 -34166.241 6279.462 304.70508 -1538.0247
40 -27892.698 3186.4294 144409.85 -182417.01 337.80177 -34097.62 6204.9214 301.08807 -1516.1201
50 -27898.215 3198.5531 144426.3 -182405.24 336.58074 -34049.909 6151.6947 298.50529 -1349.4431
60 -27900.589 3163.4592 144400.32 -182414.85 341.17705 -34110.926 6210.3369 301.35085 -1695.3697
70 -27900.487 3223.7183 144242.71 -182409.21 341.09496 -34188.493 6288.0059 305.11967 -1493.2032
80 -27901.07 3274.244 144265.07 -182417.68 344.0409 -34177.343 6276.2725 304.55032 -1273.0263
90 -27905.672 3237.6056 144288.71 -182418.22 342.15013 -34187.814 6282.1417 304.83511 -1268.0436
SHAKE stats (type/ave/delta/count) on step 100
Bond: 16 1.09 3.78281e-07 6
Bond: 18 1.09 1.12288e-07 3
Bond: 34 1.111 7.60709e-07 10
Bond: 39 1.111 2.37855e-07 5
Bond: 43 1.111 6.00872e-07 10
Bond: 44 1.111 3.75324e-07 10
Bond: 59 1.111 1.12311e-07 2
Bond: 62 1.111 2.99471e-07 2
Bond: 63 1.111 6.10589e-07 46
Bond: 64 1.111 4.50733e-07 18
Bond: 65 1.08 2.90668e-07 16
Bond: 66 1.08 1.61592e-07 4
Bond: 67 1 5.4508e-07 10
Bond: 68 1.01 0 1
Bond: 69 1.08 4.1398e-07 5
Bond: 70 1.09 0 1
Bond: 71 1.083 0 1
Bond: 72 0.9572 1.76706e-06 6530
Bond: 75 1 3.96686e-07 10
Bond: 79 0.997 7.72922e-07 17
Bond: 81 1 1.30673e-07 4
Bond: 84 1.04 1.44551e-07 9
Bond: 87 1 0 1
Bond: 95 0.96 1.03526e-07 2
Bond: 97 1.325 3.64689e-08 3
Angle: 142 104.52 0.000130126 3265
100 -27913.241 3159.0676 144299.1 -182414.94 336.48839 -34254.412 6341.1706 307.69943 -1421.2905
Loop time of 3.49837 on 4 procs for 100 steps with 10245 atoms
Performance: 4.939 ns/day, 4.859 hours/ns, 28.585 timesteps/s, 292.851 katom-step/s
99.4% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.4572 | 2.5133 | 2.5634 | 2.6 | 71.84
Bond | 0.0040264 | 0.0050718 | 0.0069286 | 1.6 | 0.14
Kspace | 0.45979 | 0.50903 | 0.56364 | 5.6 | 14.55
Neigh | 0.42029 | 0.42034 | 0.42036 | 0.0 | 12.02
Comm | 0.013207 | 0.013292 | 0.013404 | 0.1 | 0.38
Output | 0.00026525 | 0.00029549 | 0.00038249 | 0.0 | 0.01
Modify | 0.035024 | 0.035546 | 0.03621 | 0.3 | 1.02
Other | | 0.001504 | | | 0.04
Nlocal: 2561.25 ave 2599 max 2520 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Nghost: 16491.5 ave 16541 max 16442 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 1.99855e+06 ave 2.04035e+06 max 1.95468e+06 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Total # of neighbors = 7994217
Ave neighs/atom = 780.30425
Ave special neighs/atom = 2.3664226
Neighbor list builds = 10
Dangerous builds = 0
Total wall time: 0:00:03

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

@ -54,7 +54,7 @@ AngleSPICA::AngleSPICA(LAMMPS *lmp) :
AngleSPICA::~AngleSPICA() AngleSPICA::~AngleSPICA()
{ {
if (allocated) { if (allocated && !copymode) {
memory->destroy(setflag); memory->destroy(setflag);
memory->destroy(k); memory->destroy(k);
memory->destroy(theta0); memory->destroy(theta0);

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

@ -52,7 +52,7 @@ class AngleSPICA : public Angle {
void ev_tally13(int, int, int, int, double, double, double, double, double); void ev_tally13(int, int, int, int, double, double, double, double, double);
void allocate(); virtual void allocate();
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS

7
src/CG-SPICA/pair_lj_spica_coul_long.cpp
View File

@ -55,6 +55,8 @@ PairLJSPICACoulLong::PairLJSPICACoulLong(LAMMPS *lmp) :
PairLJSPICACoulLong::~PairLJSPICACoulLong() PairLJSPICACoulLong::~PairLJSPICACoulLong()
{ {
if (copymode) return;
if (allocated) { if (allocated) {
memory->destroy(setflag); memory->destroy(setflag);
memory->destroy(lj_type); memory->destroy(lj_type);
@ -356,7 +358,7 @@ void PairLJSPICACoulLong::coeff(int narg, char **arg)
void PairLJSPICACoulLong::init_style() void PairLJSPICACoulLong::init_style()
{ {
if (!atom->q_flag) error->all(FLERR, "Pair style lj/cut/coul/long requires atom attribute q"); if (!atom->q_flag) error->all(FLERR, "Pair style lj/spica/coul/long requires atom attribute q");
neighbor->add_request(this); neighbor->add_request(this);
@ -385,7 +387,8 @@ double PairLJSPICACoulLong::init_one(int i, int j)
const int ljt = lj_type[i][j]; const int ljt = lj_type[i][j];
if (ljt == LJ_NOT_SET) error->all(FLERR, "unrecognized LJ parameter flag"); if (ljt == LJ_NOT_SET)
error->all(FLERR,"unrecognized LJ parameter flag");
double cut = MAX(cut_lj[i][j], cut_coul); double cut = MAX(cut_lj[i][j], cut_coul);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];

2
src/CG-SPICA/pair_lj_spica_coul_long.h
View File

@ -64,7 +64,7 @@ class PairLJSPICACoulLong : public Pair {
double cut_lj_global; double cut_lj_global;
double g_ewald; double g_ewald;
void allocate(); virtual void allocate();
private: private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval(); template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval();

4
src/KOKKOS/Install.sh
View File

@ -55,6 +55,8 @@ action angle_harmonic_kokkos.cpp angle_harmonic.cpp
action angle_harmonic_kokkos.h angle_harmonic.h action angle_harmonic_kokkos.h angle_harmonic.h
action angle_hybrid_kokkos.cpp angle_hybrid.cpp action angle_hybrid_kokkos.cpp angle_hybrid.cpp
action angle_hybrid_kokkos.h angle_hybrid.h action angle_hybrid_kokkos.h angle_hybrid.h
action angle_spica_kokkos.cpp angle_spica.cpp
action angle_spica_kokkos.h angle_spica.h
action atom_kokkos.cpp action atom_kokkos.cpp
action atom_kokkos.h action atom_kokkos.h
action atom_map_kokkos.cpp action atom_map_kokkos.cpp
@ -350,6 +352,8 @@ action pair_lj_gromacs_coul_gromacs_kokkos.cpp pair_lj_gromacs_coul_gromacs.cpp
action pair_lj_gromacs_coul_gromacs_kokkos.h pair_lj_gromacs_coul_gromacs.h action pair_lj_gromacs_coul_gromacs_kokkos.h pair_lj_gromacs_coul_gromacs.h
action pair_lj_gromacs_kokkos.cpp pair_lj_gromacs.cpp action pair_lj_gromacs_kokkos.cpp pair_lj_gromacs.cpp
action pair_lj_gromacs_kokkos.h pair_lj_gromacs.h action pair_lj_gromacs_kokkos.h pair_lj_gromacs.h
action pair_lj_spica_coul_long_kokkos.cpp pair_lj_spica_coul_long.cpp
action pair_lj_spica_coul_long_kokkos.h pair_lj_spica_coul_long.h
action pair_lj_spica_kokkos.cpp pair_lj_spica.cpp action pair_lj_spica_kokkos.cpp pair_lj_spica.cpp
action pair_lj_spica_kokkos.h pair_lj_spica.h action pair_lj_spica_kokkos.h pair_lj_spica.h
action pair_meam_kokkos.cpp pair_meam.cpp action pair_meam_kokkos.cpp pair_meam.cpp

656
src/KOKKOS/angle_spica_kokkos.cpp Normal file
View File

@ -0,0 +1,656 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mitch Murphy (alphataubio@gmail.com)
------------------------------------------------------------------------- */
#include "angle_spica_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory_kokkos.h"
#include "neighbor_kokkos.h"
#include "respa.h"
#include "update.h"
#include "lj_spica_common.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace LJSPICAParms;
static constexpr double SMALL = 0.001;
/* ---------------------------------------------------------------------- */
template<class DeviceType>
AngleSPICAKokkos<DeviceType>::AngleSPICAKokkos(LAMMPS *lmp) : AngleSPICA(lmp)
{
atomKK = (AtomKokkos *) atom;
neighborKK = (NeighborKokkos *) neighbor;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
centroidstressflag = CENTROID_NOTAVAIL;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
AngleSPICAKokkos<DeviceType>::~AngleSPICAKokkos()
{
if (!copymode) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->destroy_kokkos(k_vatom,vatom);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void AngleSPICAKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
eflag = eflag_in;
vflag = vflag_in;
ev_init(eflag,vflag,0);
// reallocate per-atom arrays if necessary
if (eflag_atom) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"angle:eatom");
d_eatom = k_eatom.template view<DeviceType>();
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"angle:vatom");
d_vatom = k_vatom.template view<DeviceType>();
}
k_k.template sync<DeviceType>();
k_theta0.template sync<DeviceType>();
k_repscale.template sync<DeviceType>();
k_lj_type.template sync<DeviceType>();
k_lj1.template sync<DeviceType>();
k_lj2.template sync<DeviceType>();
k_lj3.template sync<DeviceType>();
k_lj4.template sync<DeviceType>();
k_rminsq.template sync<DeviceType>();
k_emin.template sync<DeviceType>();
// "It has to do with overlapping host/device in verlet_kokkos.cpp. For this reason, all topology styles (bond, angle, etc.) must set DATAMASK_READ, DATAMASK_MODIFY in the constructor and must not use atomKK->sync/modified. This is a gotcha that needed to be better documented."
// https://matsci.org/t/a-few-kokkos-development-questions/56598
//
// atomKK->sync(execution_space,datamask_read);
// if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
// else atomKK->modified(execution_space,F_MASK);
//atomKK->k_type.template sync<DeviceType>();
x = atomKK->k_x.template view<DeviceType>();
f = atomKK->k_f.template view<DeviceType>();
neighborKK->k_anglelist.template sync<DeviceType>();
anglelist = neighborKK->k_anglelist.template view<DeviceType>();
int nanglelist = neighborKK->nanglelist;
d_type = atomKK->k_type.template view<DeviceType>();
nlocal = atom->nlocal;
newton_bond = force->newton_bond;
copymode = 1;
// loop over neighbors of my atoms
EV_FLOAT ev;
if (evflag) {
if (newton_bond) {
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagAngleSPICACompute<1,1> >(0,nanglelist),*this,ev);
} else {
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagAngleSPICACompute<0,1> >(0,nanglelist),*this,ev);
}
} else {
if (newton_bond) {
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagAngleSPICACompute<1,0> >(0,nanglelist),*this);
} else {
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagAngleSPICACompute<0,0> >(0,nanglelist),*this);
}
}
if (eflag_global) energy += ev.evdwl;
if (vflag_global) {
virial[0] += ev.v[0];
virial[1] += ev.v[1];
virial[2] += ev.v[2];
virial[3] += ev.v[3];
virial[4] += ev.v[4];
virial[5] += ev.v[5];
}
if (eflag_atom) {
k_eatom.template modify<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
}
template<class DeviceType>
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void AngleSPICAKokkos<DeviceType>::operator()(TagAngleSPICACompute<NEWTON_BOND,EVFLAG>, const int &n, EV_FLOAT& ev) const {
// The f array is atomic
Kokkos::View<F_FLOAT*[3], typename DAT::t_f_array::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<Kokkos::Atomic|Kokkos::Unmanaged> > a_f = f;
const int i1 = anglelist(n,0);
const int i2 = anglelist(n,1);
const int i3 = anglelist(n,2);
const int type = anglelist(n,3);
// 1st bond
const F_FLOAT delx1 = x(i1,0) - x(i2,0);
const F_FLOAT dely1 = x(i1,1) - x(i2,1);
const F_FLOAT delz1 = x(i1,2) - x(i2,2);
const F_FLOAT rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
const F_FLOAT r1 = sqrt(rsq1);
// 2nd bond
const F_FLOAT delx2 = x(i3,0) - x(i2,0);
const F_FLOAT dely2 = x(i3,1) - x(i2,1);
const F_FLOAT delz2 = x(i3,2) - x(i2,2);
const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
const F_FLOAT r2 = sqrt(rsq2);
// angle (cos and sin)
F_FLOAT c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
F_FLOAT s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL;
s = 1.0/s;
// 1-3 LJ interaction.
// we only want to use the repulsive part,
// and it can be scaled (or off).
// so this has to be done here and not in the
// general non-bonded code.
F_FLOAT f13, e13, delx3, dely3, delz3;
f13 = e13 = delx3 = dely3 = delz3 = 0.0;
if (repflag) {
delx3 = x(i1,0) - x(i3,0);
dely3 = x(i1,1) - x(i3,1);
delz3 = x(i1,2) - x(i3,2);
const F_FLOAT rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3;
const int type1 = d_type[i1];
const int type3 = d_type[i3];
f13=0.0;
e13=0.0;
if (rsq3 < d_rminsq(type1,type3)) {
const int ljt = d_lj_type(type1,type3);
const double r2inv = 1.0/rsq3;
if (ljt == LJ12_4) {
const double r4inv=r2inv*r2inv;
f13 = r4inv*(d_lj1(type1,type3)*r4inv*r4inv - d_lj2(type1,type3));
if (eflag) e13 = r4inv*(d_lj3(type1,type3)*r4inv*r4inv - d_lj4(type1,type3));
} else if (ljt == LJ9_6) {
const double r3inv = r2inv*sqrt(r2inv);
const double r6inv = r3inv*r3inv;
f13 = r6inv*(d_lj1(type1,type3)*r3inv - d_lj2(type1,type3));
if (eflag) e13 = r6inv*(d_lj3(type1,type3)*r3inv - d_lj4(type1,type3));
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
f13 = r6inv*(d_lj1(type1,type3)*r6inv - d_lj2(type1,type3));
if (eflag) e13 = r6inv*(d_lj3(type1,type3)*r6inv - d_lj4(type1,type3));
} else if (ljt == LJ12_5) {
const double r5inv = r2inv*r2inv*sqrt(r2inv);
const double r7inv = r5inv*r2inv;
f13 = r5inv*(d_lj1(type1,type3)*r7inv - d_lj2(type1,type3));
if (eflag) e13 = r5inv*(d_lj3(type1,type3)*r7inv - d_lj4(type1,type3));
}
// make sure energy is 0.0 at the cutoff.
if (eflag) e13 -= d_emin(type1,type3);
f13 *= r2inv;
}
}
// force & energy
const F_FLOAT dtheta = acos(c) - d_theta0[type];
const F_FLOAT tk = d_k[type] * dtheta;
F_FLOAT eangle = 0.0;
if (eflag) eangle = tk*dtheta;
const F_FLOAT a = -2.0 * tk * s;
const F_FLOAT a11 = a*c / rsq1;
const F_FLOAT a12 = -a / (r1*r2);
const F_FLOAT a22 = a*c / rsq2;
F_FLOAT f1[3],f3[3];
f1[0] = a11*delx1 + a12*delx2;
f1[1] = a11*dely1 + a12*dely2;
f1[2] = a11*delz1 + a12*delz2;
f3[0] = a22*delx2 + a12*delx1;
f3[1] = a22*dely2 + a12*dely1;
f3[2] = a22*delz2 + a12*delz1;
// apply force to each of 3 atoms
if (NEWTON_BOND || i1 < nlocal) {
a_f(i1,0) += f1[0] + f13*delx3;
a_f(i1,1) += f1[1] + f13*dely3;
a_f(i1,2) += f1[2] + f13*delz3;
}
if (NEWTON_BOND || i2 < nlocal) {
a_f(i2,0) -= f1[0] + f3[0];
a_f(i2,1) -= f1[1] + f3[1];
a_f(i2,2) -= f1[2] + f3[2];
}
if (NEWTON_BOND || i3 < nlocal) {
a_f(i3,0) += f3[0] - f13*delx3;
a_f(i3,1) += f3[1] - f13*dely3;
a_f(i3,2) += f3[2] - f13*delz3;
}
if (EVFLAG) {
ev_tally(ev,i1,i2,i3,eangle,f1,f3,delx1,dely1,delz1,delx2,dely2,delz2);
if (repflag)
ev_tally13(ev,i1,i3,e13,f13,delx3,dely3,delz3);
}
}
template<class DeviceType>
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void AngleSPICAKokkos<DeviceType>::operator()(TagAngleSPICACompute<NEWTON_BOND,EVFLAG>, const int &n) const {
EV_FLOAT ev;
this->template operator()<NEWTON_BOND,EVFLAG>(TagAngleSPICACompute<NEWTON_BOND,EVFLAG>(), n, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void AngleSPICAKokkos<DeviceType>::allocate()
{
AngleSPICA::allocate();
int nangletypes = atom->nangletypes;
k_k = typename ArrayTypes<DeviceType>::tdual_ffloat_1d("AngleSPICA::k",nangletypes+1);
k_theta0 = typename ArrayTypes<DeviceType>::tdual_ffloat_1d("AngleSPICA::theta0",nangletypes+1);
k_repscale = typename ArrayTypes<DeviceType>::tdual_ffloat_1d("AngleSPICA::repscale",nangletypes+1);
k_setflag = typename ArrayTypes<DeviceType>::tdual_int_1d("AngleSPICA::setflag",nangletypes+1);
d_k = k_k.template view<DeviceType>();
d_theta0 = k_theta0.template view<DeviceType>();
d_repscale = k_repscale.template view<DeviceType>();
d_setflag = k_setflag.template view<DeviceType>();
int ntypes = atom->ntypes;
k_lj_type = typename ArrayTypes<DeviceType>::tdual_int_2d("AngleSPICA::lj_type",ntypes+1,ntypes+1);
k_lj1 = typename ArrayTypes<DeviceType>::tdual_ffloat_2d("AngleSPICA::lj1",ntypes+1,ntypes+1);
k_lj2 = typename ArrayTypes<DeviceType>::tdual_ffloat_2d("AngleSPICA::lj2",ntypes+1,ntypes+1);
k_lj3 = typename ArrayTypes<DeviceType>::tdual_ffloat_2d("AngleSPICA::lj3",ntypes+1,ntypes+1);
k_lj4 = typename ArrayTypes<DeviceType>::tdual_ffloat_2d("AngleSPICA::lj4",ntypes+1,ntypes+1);
k_rminsq = typename ArrayTypes<DeviceType>::tdual_ffloat_2d("AngleSPICA::rminsq",ntypes+1,ntypes+1);
k_emin = typename ArrayTypes<DeviceType>::tdual_ffloat_2d("AngleSPICA::emin",ntypes+1,ntypes+1);
d_lj_type = k_lj_type.template view<DeviceType>();
d_lj1 = k_lj1.template view<DeviceType>();
d_lj2 = k_lj2.template view<DeviceType>();
d_lj3 = k_lj3.template view<DeviceType>();
d_lj4 = k_lj4.template view<DeviceType>();
d_rminsq = k_rminsq.template view<DeviceType>();
d_emin = k_emin.template view<DeviceType>();
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
template<class DeviceType>
void AngleSPICAKokkos<DeviceType>::init_style()
{
AngleSPICA::init_style();
// error if rRESPA with inner levels
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
int respa = 0;
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
if (respa)
error->all(FLERR,"Cannot use Kokkos pair style with rRESPA inner/middle");
}
int ntypes = atom->ntypes;
for (int i = 1; i <= ntypes; i++) {
for (int j = 1; j <= ntypes; j++) {
k_lj_type.h_view(i,j) = lj_type[i][j];
k_lj1.h_view(i,j) = lj1[i][j];
k_lj2.h_view(i,j) = lj2[i][j];
k_lj3.h_view(i,j) = lj3[i][j];
k_lj4.h_view(i,j) = lj4[i][j];
k_rminsq.h_view(i,j) = rminsq[i][j];
k_emin.h_view(i,j) = emin[i][j];
}
}
k_lj_type.template modify<LMPHostType>();
k_lj1.template modify<LMPHostType>();
k_lj2.template modify<LMPHostType>();
k_lj3.template modify<LMPHostType>();
k_lj4.template modify<LMPHostType>();
k_rminsq.template modify<LMPHostType>();
k_emin.template modify<LMPHostType>();
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
template<class DeviceType>
void AngleSPICAKokkos<DeviceType>::coeff(int narg, char **arg)
{
AngleSPICA::coeff(narg, arg);
int n = atom->nangletypes;
for (int i = 1; i <= n; i++) {
k_k.h_view[i] = k[i];
k_theta0.h_view[i] = theta0[i];
k_repscale.h_view[i] = repscale[i];
k_setflag.h_view[i] = setflag[i];
}
k_k.template modify<LMPHostType>();
k_theta0.template modify<LMPHostType>();
k_repscale.template modify<LMPHostType>();
k_setflag.template modify<LMPHostType>();
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
template<class DeviceType>
void AngleSPICAKokkos<DeviceType>::read_restart(FILE *fp)
{
AngleSPICA::read_restart(fp);
int n = atom->nangletypes;
for (int i = 1; i <= n; i++) {
k_k.h_view[i] = k[i];
k_theta0.h_view[i] = theta0[i];
k_repscale.h_view[i] = repscale[i];
k_setflag.h_view[i] = setflag[i];
}
k_k.template modify<LMPHostType>();
k_theta0.template modify<LMPHostType>();
k_repscale.template modify<LMPHostType>();
k_setflag.template modify<LMPHostType>();
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
------------------------------------------------------------------------- */
template<class DeviceType>
//template<int NEWTON_BOND>
KOKKOS_INLINE_FUNCTION
void AngleSPICAKokkos<DeviceType>::ev_tally(EV_FLOAT &ev, const int i, const int j, const int k,
F_FLOAT &eangle, F_FLOAT *f1, F_FLOAT *f3,
const F_FLOAT &delx1, const F_FLOAT &dely1, const F_FLOAT &delz1,
const F_FLOAT &delx2, const F_FLOAT &dely2, const F_FLOAT &delz2) const
{
E_FLOAT eanglethird;
F_FLOAT v[6];
// The eatom and vatom arrays are atomic
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<Kokkos::Atomic|Kokkos::Unmanaged> > v_eatom = k_eatom.template view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<Kokkos::Atomic|Kokkos::Unmanaged> > v_vatom = k_vatom.template view<DeviceType>();
if (eflag_either) {
if (eflag_global) {
if (newton_bond) ev.evdwl += eangle;
else {
eanglethird = THIRD*eangle;
if (i < nlocal) ev.evdwl += eanglethird;
if (j < nlocal) ev.evdwl += eanglethird;
if (k < nlocal) ev.evdwl += eanglethird;
}
}
if (eflag_atom) {
eanglethird = THIRD*eangle;
if (newton_bond || i < nlocal) v_eatom[i] += eanglethird;
if (newton_bond || j < nlocal) v_eatom[j] += eanglethird;
if (newton_bond || k < nlocal) v_eatom[k] += eanglethird;
}
}
if (vflag_either) {
v[0] = delx1*f1[0] + delx2*f3[0];
v[1] = dely1*f1[1] + dely2*f3[1];
v[2] = delz1*f1[2] + delz2*f3[2];
v[3] = delx1*f1[1] + delx2*f3[1];
v[4] = delx1*f1[2] + delx2*f3[2];
v[5] = dely1*f1[2] + dely2*f3[2];
if (vflag_global) {
if (newton_bond) {
ev.v[0] += v[0];
ev.v[1] += v[1];
ev.v[2] += v[2];
ev.v[3] += v[3];
ev.v[4] += v[4];
ev.v[5] += v[5];
} else {
if (i < nlocal) {
ev.v[0] += THIRD*v[0];
ev.v[1] += THIRD*v[1];
ev.v[2] += THIRD*v[2];
ev.v[3] += THIRD*v[3];
ev.v[4] += THIRD*v[4];
ev.v[5] += THIRD*v[5];
}
if (j < nlocal) {
ev.v[0] += THIRD*v[0];
ev.v[1] += THIRD*v[1];
ev.v[2] += THIRD*v[2];
ev.v[3] += THIRD*v[3];
ev.v[4] += THIRD*v[4];
ev.v[5] += THIRD*v[5];
}
if (k < nlocal) {
ev.v[0] += THIRD*v[0];
ev.v[1] += THIRD*v[1];
ev.v[2] += THIRD*v[2];
ev.v[3] += THIRD*v[3];
ev.v[4] += THIRD*v[4];
ev.v[5] += THIRD*v[5];
}
}
}
if (vflag_atom) {
if (newton_bond || i < nlocal) {
v_vatom(i,0) += THIRD*v[0];
v_vatom(i,1) += THIRD*v[1];
v_vatom(i,2) += THIRD*v[2];
v_vatom(i,3) += THIRD*v[3];
v_vatom(i,4) += THIRD*v[4];
v_vatom(i,5) += THIRD*v[5];
}
if (newton_bond || j < nlocal) {
v_vatom(j,0) += THIRD*v[0];
v_vatom(j,1) += THIRD*v[1];
v_vatom(j,2) += THIRD*v[2];
v_vatom(j,3) += THIRD*v[3];
v_vatom(j,4) += THIRD*v[4];
v_vatom(j,5) += THIRD*v[5];
}
if (newton_bond || k < nlocal) {
v_vatom(k,0) += THIRD*v[0];
v_vatom(k,1) += THIRD*v[1];
v_vatom(k,2) += THIRD*v[2];
v_vatom(k,3) += THIRD*v[3];
v_vatom(k,4) += THIRD*v[4];
v_vatom(k,5) += THIRD*v[5];
}
}
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void AngleSPICAKokkos<DeviceType>::ev_tally13(EV_FLOAT &ev, const int i, const int j,
const F_FLOAT &evdwl, const F_FLOAT &fpair,
const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const
{
double v[6];
// The eatom and vatom arrays are atomic
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<Kokkos::Atomic|Kokkos::Unmanaged> > v_eatom = k_eatom.template view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,typename KKDevice<DeviceType>::value,Kokkos::MemoryTraits<Kokkos::Atomic|Kokkos::Unmanaged> > v_vatom = k_vatom.template view<DeviceType>();
if (eflag_either) {
if (eflag_global) {
if (newton_bond) {
ev.evdwl += evdwl;
} else {
if (i < nlocal)
ev.evdwl += 0.5*evdwl;
if (j < nlocal)
ev.evdwl += 0.5*evdwl;
}
}
if (eflag_atom) {
if (newton_bond || i < nlocal) v_eatom[i] += 0.5*evdwl;
if (newton_bond || j < nlocal) v_eatom[j] += 0.5*evdwl;
}
}
if (vflag_either) {
v[0] = delx*delx*fpair;
v[1] = dely*dely*fpair;
v[2] = delz*delz*fpair;
v[3] = delx*dely*fpair;
v[4] = delx*delz*fpair;
v[5] = dely*delz*fpair;
if (vflag_global) {
if (newton_bond) {
ev.v[0] += v[0];
ev.v[1] += v[1];
ev.v[2] += v[2];
ev.v[3] += v[3];
ev.v[4] += v[4];
ev.v[5] += v[5];
} else {
if (i < nlocal) {
ev.v[0] += 0.5*v[0];
ev.v[1] += 0.5*v[1];
ev.v[2] += 0.5*v[2];
ev.v[3] += 0.5*v[3];
ev.v[4] += 0.5*v[4];
ev.v[5] += 0.5*v[5];
}
if (j < nlocal) {
ev.v[0] += 0.5*v[0];
ev.v[1] += 0.5*v[1];
ev.v[2] += 0.5*v[2];
ev.v[3] += 0.5*v[3];
ev.v[4] += 0.5*v[4];
ev.v[5] += 0.5*v[5];
}
}
}
if (vflag_atom) {
if (newton_bond || i < nlocal) {
v_vatom(i,0) += 0.5*v[0];
v_vatom(i,1) += 0.5*v[1];
v_vatom(i,2) += 0.5*v[2];
v_vatom(i,3) += 0.5*v[3];
v_vatom(i,4) += 0.5*v[4];
v_vatom(i,5) += 0.5*v[5];
}
if (newton_bond || j < nlocal) {
v_vatom(j,0) += 0.5*v[0];
v_vatom(j,1) += 0.5*v[1];
v_vatom(j,2) += 0.5*v[2];
v_vatom(j,3) += 0.5*v[3];
v_vatom(j,4) += 0.5*v[4];
v_vatom(j,5) += 0.5*v[5];
}
}
}
}
/* ---------------------------------------------------------------------- */
namespace LAMMPS_NS {
template class AngleSPICAKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class AngleSPICAKokkos<LMPHostType>;
#endif
}

106
src/KOKKOS/angle_spica_kokkos.h Normal file
View File

@ -0,0 +1,106 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef ANGLE_CLASS
// clang-format off
AngleStyle(spica/kk,AngleSPICAKokkos<LMPDeviceType>);
AngleStyle(spica/kk/device,AngleSPICAKokkos<LMPDeviceType>);
AngleStyle(spica/kk/host,AngleSPICAKokkos<LMPHostType>);
AngleStyle(sdk/kk,AngleSPICAKokkos<LMPDeviceType>);
AngleStyle(sdk/kk/device,AngleSPICAKokkos<LMPDeviceType>);
AngleStyle(sdk/kk/host,AngleSPICAKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_ANGLE_SPICA_KOKKOS_H
#define LMP_ANGLE_SPICA_KOKKOS_H
#include "angle_spica.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
template<int NEWTON_BOND, int EVFLAG>
struct TagAngleSPICACompute{};
template<class DeviceType>
class AngleSPICAKokkos : public AngleSPICA {
public:
typedef DeviceType device_type;
typedef EV_FLOAT value_type;
AngleSPICAKokkos(class LAMMPS *);
~AngleSPICAKokkos() override;
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
void read_restart(FILE *) override;
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagAngleSPICACompute<NEWTON_BOND,EVFLAG>, const int&, EV_FLOAT&) const;
template<int NEWTON_BOND, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagAngleSPICACompute<NEWTON_BOND,EVFLAG>, const int&) const;
//template<int NEWTON_BOND>
KOKKOS_INLINE_FUNCTION
void ev_tally(EV_FLOAT &ev, const int i, const int j, const int k,
F_FLOAT &eangle, F_FLOAT *f1, F_FLOAT *f3,
const F_FLOAT &delx1, const F_FLOAT &dely1, const F_FLOAT &delz1,
const F_FLOAT &delx2, const F_FLOAT &dely2, const F_FLOAT &delz2) const;
KOKKOS_INLINE_FUNCTION
void ev_tally13(EV_FLOAT &ev, const int i, const int j,
const F_FLOAT &evdwl, const F_FLOAT &fpair,
const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
protected:
class NeighborKokkos *neighborKK;
typename ArrayTypes<DeviceType>::t_x_array_randomread x;
typename ArrayTypes<DeviceType>::t_f_array f;
typename ArrayTypes<DeviceType>::t_int_2d anglelist;
typename ArrayTypes<DeviceType>::tdual_efloat_1d k_eatom;
typename ArrayTypes<DeviceType>::tdual_virial_array k_vatom;
typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
int nlocal,newton_bond;
int eflag,vflag;
typename ArrayTypes<DeviceType>::tdual_int_1d k_setflag;
typename ArrayTypes<DeviceType>::t_int_1d d_setflag, d_type;
typename ArrayTypes<DeviceType>::tdual_ffloat_1d k_k, k_theta0, k_repscale;
typename ArrayTypes<DeviceType>::t_ffloat_1d d_k, d_theta0, d_repscale;
typename ArrayTypes<DeviceType>::tdual_int_2d k_lj_type;
typename ArrayTypes<DeviceType>::t_int_2d d_lj_type;
typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_lj1, k_lj2, k_lj3, k_lj4, k_rminsq, k_emin;
typename ArrayTypes<DeviceType>::t_ffloat_2d d_lj1, d_lj2, d_lj3, d_lj4, d_rminsq, d_emin;
void allocate() override;
};
}
#endif
#endif

16
src/KOKKOS/atom_kokkos.cpp
View File

@ -276,22 +276,6 @@ void AtomKokkos::sort_device()
if (domain->triclinic) domain->x2lamda(nlocal); if (domain->triclinic) domain->x2lamda(nlocal);
} }
/* ----------------------------------------------------------------------
reallocate memory to the pointer selected by the mask
------------------------------------------------------------------------- */
void AtomKokkos::grow(unsigned int mask)
{
if (mask & SPECIAL_MASK) {
memoryKK->destroy_kokkos(k_special, special);
sync(Device, mask);
modified(Device, mask);
memoryKK->grow_kokkos(k_special, special, nmax, maxspecial, "atom:special");
avec->grow_pointers();
sync(Host, mask);
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add a custom variable with name of type flag = 0/1 for int/double add a custom variable with name of type flag = 0/1 for int/double
assumes name does not already exist assumes name does not already exist

1
src/KOKKOS/atom_kokkos.h
View File

@ -165,7 +165,6 @@ class AtomKokkos : public Atom {
void modified(const ExecutionSpace space, unsigned int mask); void modified(const ExecutionSpace space, unsigned int mask);
void sync_overlapping_device(const ExecutionSpace space, unsigned int mask); void sync_overlapping_device(const ExecutionSpace space, unsigned int mask);
void sort() override; void sort() override;
virtual void grow(unsigned int mask);
int add_custom(const char *, int, int, int border = 0) override; int add_custom(const char *, int, int, int border = 0) override;
void remove_custom(int, int, int) override; void remove_custom(int, int, int) override;
void deallocate_topology() override; void deallocate_topology() override;

3
src/KOKKOS/pair_lj_charmmfsw_coul_long_kokkos.cpp
View File

@ -13,7 +13,7 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Contributing author: Mitch Murphy (alphataubio) Contributing author: Mitch Murphy (alphataubio@gmail.com)
Based on serial kspace lj-fsw sections (force-switched) provided by Based on serial kspace lj-fsw sections (force-switched) provided by
Robert Meissner and Lucio Colombi Ciacchi of Bremen University, Germany, Robert Meissner and Lucio Colombi Ciacchi of Bremen University, Germany,
@ -463,7 +463,6 @@ double PairLJCharmmfswCoulLongKokkos<DeviceType>::init_one(int i, int j)
k_params.h_view(i,j).lj2 = lj2[i][j]; k_params.h_view(i,j).lj2 = lj2[i][j];
k_params.h_view(i,j).lj3 = lj3[i][j]; k_params.h_view(i,j).lj3 = lj3[i][j];
k_params.h_view(i,j).lj4 = lj4[i][j]; k_params.h_view(i,j).lj4 = lj4[i][j];
//k_params.h_view(i,j).offset = offset[i][j];
k_params.h_view(i,j).cut_ljsq = cut_ljsq; k_params.h_view(i,j).cut_ljsq = cut_ljsq;
k_params.h_view(i,j).cut_coulsq = cut_coulsq; k_params.h_view(i,j).cut_coulsq = cut_coulsq;

27
src/KOKKOS/pair_lj_charmmfsw_coul_long_kokkos.h
View File

@ -48,30 +48,25 @@ class PairLJCharmmfswCoulLongKokkos : public PairLJCharmmfswCoulLong {
protected: protected:
template<bool STACKPARAMS, class Specialisation> template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fpair(const F_FLOAT& rsq, const int& i, const int&j, F_FLOAT compute_fpair(const F_FLOAT& rsq, const int& i, const int& j, const int& itype, const int& jtype) const;
const int& itype, const int& jtype) const;
template<bool STACKPARAMS, class Specialisation> template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fcoul(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, F_FLOAT compute_evdwl(const F_FLOAT& rsq, const int& i, const int& j, const int& itype, const int& jtype) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fcoul(const F_FLOAT& rsq, const int& i, const int& j, const int& itype,
const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const; const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
template<bool STACKPARAMS, class Specialisation> template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
F_FLOAT compute_evdwl(const F_FLOAT& rsq, const int& i, const int&j, F_FLOAT compute_ecoul(const F_FLOAT& rsq, const int& i, const int& j, const int& itype,
const int& itype, const int& jtype) const; const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_ecoul(const F_FLOAT& rsq, const int& i, const int&j,
const int& itype, const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params; Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_lj_coul**, typename Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
Kokkos::LayoutRight,DeviceType>::t_dev_const_um params; params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1]; // hardwired to space for 12 atom types
// hardwired to space for 12 atom types
params_lj_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1]; F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1]; F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1]; F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
@ -100,8 +95,8 @@ class PairLJCharmmfswCoulLongKokkos : public PairLJCharmmfswCoulLong {
int neighflag; int neighflag;
int nlocal,nall,eflag,vflag; int nlocal,nall,eflag,vflag;
double special_coul[4];
double special_lj[4]; double special_lj[4];
double special_coul[4];
double qqrd2e; double qqrd2e;
void allocate() override; void allocate() override;

3
src/KOKKOS/pair_lj_cut_coul_long_kokkos.cpp
View File

@ -76,7 +76,6 @@ void PairLJCutCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
ev_init(eflag,vflag,0); ev_init(eflag,vflag,0);
// reallocate per-atom arrays if necessary // reallocate per-atom arrays if necessary
if (eflag_atom) { if (eflag_atom) {
@ -125,11 +124,11 @@ void PairLJCutCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
ev = pair_compute<PairLJCutCoulLongKokkos<DeviceType>,CoulLongTable<0> > ev = pair_compute<PairLJCutCoulLongKokkos<DeviceType>,CoulLongTable<0> >
(this,(NeighListKokkos<DeviceType>*)list); (this,(NeighListKokkos<DeviceType>*)list);
if (eflag) { if (eflag) {
eng_vdwl += ev.evdwl; eng_vdwl += ev.evdwl;
eng_coul += ev.ecoul; eng_coul += ev.ecoul;
} }
if (vflag_global) { if (vflag_global) {
virial[0] += ev.v[0]; virial[0] += ev.v[0];
virial[1] += ev.v[1]; virial[1] += ev.v[1];

503
src/KOKKOS/pair_lj_spica_coul_long_kokkos.cpp Normal file
View File

@ -0,0 +1,503 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mitch Murphy (alphataubio@gmail.com)
------------------------------------------------------------------------- */
#include "pair_lj_spica_coul_long_kokkos.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "error.h"
#include "ewald_const.h"
#include "force.h"
#include "kokkos.h"
#include "memory_kokkos.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
#include "lj_spica_common.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace LJSPICAParms;
using namespace EwaldConst;
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairLJSPICACoulLongKokkos<DeviceType>::PairLJSPICACoulLongKokkos(LAMMPS *lmp) : PairLJSPICACoulLong(lmp)
{
respa_enable = 0;
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | TYPE_MASK | Q_MASK | ENERGY_MASK | VIRIAL_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairLJSPICACoulLongKokkos<DeviceType>::~PairLJSPICACoulLongKokkos()
{
if (copymode) return;
if (allocated) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->destroy_kokkos(k_cutsq,cutsq);
memoryKK->destroy_kokkos(k_cut_ljsq,cut_ljsq);
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairLJSPICACoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
eflag = eflag_in;
vflag = vflag_in;
if (neighflag == FULL) no_virial_fdotr_compute = 1;
ev_init(eflag,vflag,0);
// reallocate per-atom arrays if necessary
if (eflag_atom) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.view<DeviceType>();
}
if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
d_vatom = k_vatom.view<DeviceType>();
}
atomKK->sync(execution_space,datamask_read);
k_cutsq.template sync<DeviceType>();
k_cut_ljsq.template sync<DeviceType>();
k_params.template sync<DeviceType>();
if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
else atomKK->modified(execution_space,F_MASK);
x = atomKK->k_x.view<DeviceType>();
c_x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
q = atomKK->k_q.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
special_lj[0] = force->special_lj[0];
special_lj[1] = force->special_lj[1];
special_lj[2] = force->special_lj[2];
special_lj[3] = force->special_lj[3];
special_coul[0] = force->special_coul[0];
special_coul[1] = force->special_coul[1];
special_coul[2] = force->special_coul[2];
special_coul[3] = force->special_coul[3];
qqrd2e = force->qqrd2e;
newton_pair = force->newton_pair;
// loop over neighbors of my atoms
copymode = 1;
EV_FLOAT ev;
if (ncoultablebits)
ev = pair_compute<PairLJSPICACoulLongKokkos<DeviceType>,CoulLongTable<1> >
(this,(NeighListKokkos<DeviceType>*)list);
else
ev = pair_compute<PairLJSPICACoulLongKokkos<DeviceType>,CoulLongTable<0> >
(this,(NeighListKokkos<DeviceType>*)list);
if (eflag) {
eng_vdwl += ev.evdwl;
eng_coul += ev.ecoul;
}
if (vflag_global) {
virial[0] += ev.v[0];
virial[1] += ev.v[1];
virial[2] += ev.v[2];
virial[3] += ev.v[3];
virial[4] += ev.v[4];
virial[5] += ev.v[5];
}
if (eflag_atom) {
k_eatom.template modify<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
if (vflag_fdotr) pair_virial_fdotr_compute(this);
copymode = 0;
}
/* ----------------------------------------------------------------------
compute pair force between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJSPICACoulLongKokkos<DeviceType>::
compute_fpair(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, const int& jtype) const {
(void) i;
(void) j;
const F_FLOAT r2inv = 1.0/rsq;
const int ljt = (STACKPARAMS?m_params[itype][jtype].lj_type:params(itype,jtype).lj_type);
const F_FLOAT lj_1 = (STACKPARAMS?m_params[itype][jtype].lj1:params(itype,jtype).lj1);
const F_FLOAT lj_2 = (STACKPARAMS?m_params[itype][jtype].lj2:params(itype,jtype).lj2);
const F_FLOAT r4inv=r2inv*r2inv;
const F_FLOAT r6inv=r2inv*r4inv;
const F_FLOAT a = ljt==LJ12_4?r4inv:(ljt==LJ12_5?r4inv*sqrt(r2inv):r6inv);
const F_FLOAT b = ljt==LJ12_4?r4inv:(ljt==LJ9_6?1.0/sqrt(r2inv):(ljt==LJ12_5?r2inv*sqrt(r2inv):r2inv));
return a* ( lj_1*r6inv*b - lj_2 * r2inv);
}
/* ----------------------------------------------------------------------
compute pair potential energy between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJSPICACoulLongKokkos<DeviceType>::
compute_evdwl(const F_FLOAT& rsq, const int& i, const int&j, const int& itype, const int& jtype) const {
(void) i;
(void) j;
const F_FLOAT r2inv = 1.0/rsq;
const int ljt = (STACKPARAMS?m_params[itype][jtype].lj_type:params(itype,jtype).lj_type);
const F_FLOAT lj_3 = (STACKPARAMS?m_params[itype][jtype].lj3:params(itype,jtype).lj3);
const F_FLOAT lj_4 = (STACKPARAMS?m_params[itype][jtype].lj4:params(itype,jtype).lj4);
const F_FLOAT offset = (STACKPARAMS?m_params[itype][jtype].offset:params(itype,jtype).offset);
if (ljt == LJ12_4) {
const F_FLOAT r4inv=r2inv*r2inv;
return r4inv*(lj_3*r4inv*r4inv - lj_4) - offset;
} else if (ljt == LJ9_6) {
const F_FLOAT r3inv = r2inv*sqrt(r2inv);
const F_FLOAT r6inv = r3inv*r3inv;
return r6inv*(lj_3*r3inv - lj_4) - offset;
} else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv;
return r6inv*(lj_3*r6inv - lj_4) - offset;
} else if (ljt == LJ12_5) {
const F_FLOAT r5inv = r2inv*r2inv*sqrt(r2inv);
const F_FLOAT r7inv = r5inv*r2inv;
return r5inv*(lj_3*r7inv - lj_4) - offset;
} else
return 0.0;
}
/* ----------------------------------------------------------------------
compute coulomb pair force between atoms i and j
------------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJSPICACoulLongKokkos<DeviceType>::
compute_fcoul(const F_FLOAT& rsq, const int& /*i*/, const int&j,
const int& /*itype*/, const int& /*jtype*/,
const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const {
if (Specialisation::DoTable && rsq > tabinnersq) {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
const F_FLOAT table = d_ftable[itable] + fraction*d_dftable[itable];
F_FLOAT forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
const F_FLOAT table = d_ctable[itable] + fraction*d_dctable[itable];
const F_FLOAT prefactor = qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
return forcecoul/rsq;
} else {
const F_FLOAT r = sqrt(rsq);
const F_FLOAT grij = g_ewald * r;
const F_FLOAT expm2 = exp(-grij*grij);
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
const F_FLOAT rinv = 1.0/r;
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]*rinv;
F_FLOAT forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
return forcecoul*rinv*rinv;
}
}
/* ----------------------------------------------------------------------
compute coulomb pair potential energy between atoms i and j
------------------------------------------------------------------------- */
template<class DeviceType>
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT PairLJSPICACoulLongKokkos<DeviceType>::
compute_ecoul(const F_FLOAT& rsq, const int& /*i*/, const int&j,
const int& /*itype*/, const int& /*jtype*/, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const {
if (Specialisation::DoTable && rsq > tabinnersq) {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
const int itable = (rsq_lookup.i & ncoulmask) >> ncoulshiftbits;
const F_FLOAT fraction = (rsq_lookup.f - d_rtable[itable]) * d_drtable[itable];
const F_FLOAT table = d_etable[itable] + fraction*d_detable[itable];
F_FLOAT ecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
const F_FLOAT table = d_ctable[itable] + fraction*d_dctable[itable];
const F_FLOAT prefactor = qtmp*q[j] * table;
ecoul -= (1.0-factor_coul)*prefactor;
}
return ecoul;
} else {
const F_FLOAT r = sqrt(rsq);
const F_FLOAT grij = g_ewald * r;
const F_FLOAT expm2 = exp(-grij*grij);
const F_FLOAT t = 1.0 / (1.0 + EWALD_P*grij);
const F_FLOAT erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
const F_FLOAT prefactor = qqrd2e * qtmp*q[j]/r;
F_FLOAT ecoul = prefactor * erfc;
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
return ecoul;
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
template<class DeviceType>
void PairLJSPICACoulLongKokkos<DeviceType>::allocate()
{
PairLJSPICACoulLong::allocate();
int n = atom->ntypes;
memory->destroy(cutsq);
memoryKK->create_kokkos(k_cutsq,cutsq,n+1,n+1,"pair:cutsq");
d_cutsq = k_cutsq.template view<DeviceType>();
memory->destroy(cut_ljsq);
memoryKK->create_kokkos(k_cut_ljsq,cut_ljsq,n+1,n+1,"pair:cut_ljsq");
d_cut_ljsq = k_cut_ljsq.template view<DeviceType>();
d_cut_coulsq = typename AT::t_ffloat_2d("pair:cut_coulsq",n+1,n+1);
k_params = Kokkos::DualView<params_lj_spica_coul**,Kokkos::LayoutRight,DeviceType>("PairLJSPICACoulLong::params",n+1,n+1);
params = k_params.template view<DeviceType>();
}
/* ----------------------------------------------------------------------
init tables
------------------------------------------------------------------------- */
template<class DeviceType>
void PairLJSPICACoulLongKokkos<DeviceType>::init_tables(double cut_coul, double *cut_respa)
{
Pair::init_tables(cut_coul,cut_respa);
typedef typename ArrayTypes<DeviceType>::t_ffloat_1d table_type;
typedef typename ArrayTypes<LMPHostType>::t_ffloat_1d host_table_type;
int ntable = 1;
for (int i = 0; i < ncoultablebits; i++) ntable *= 2;
// Copy rtable and drtable
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = rtable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_rtable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = drtable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_drtable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
// Copy ftable and dftable
for (int i = 0; i < ntable; i++) {
h_table(i) = ftable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_ftable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = dftable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_dftable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
// Copy ctable and dctable
for (int i = 0; i < ntable; i++) {
h_table(i) = ctable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_ctable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = dctable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_dctable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
// Copy etable and detable
for (int i = 0; i < ntable; i++) {
h_table(i) = etable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_etable = d_table;
}
{
host_table_type h_table("HostTable",ntable);
table_type d_table("DeviceTable",ntable);
for (int i = 0; i < ntable; i++) {
h_table(i) = detable[i];
}
Kokkos::deep_copy(d_table,h_table);
d_detable = d_table;
}
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
template<class DeviceType>
void PairLJSPICACoulLongKokkos<DeviceType>::init_style()
{
PairLJSPICACoulLong::init_style();
Kokkos::deep_copy(d_cut_coulsq,cut_coulsq);
// error if rRESPA with inner levels
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
int respa = 0;
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
if (respa)
error->all(FLERR,"Cannot use Kokkos pair style with rRESPA inner/middle");
}
// adjust neighbor list request for KOKKOS
neighflag = lmp->kokkos->neighflag;
auto request = neighbor->find_request(this);
request->set_kokkos_host(std::is_same_v<DeviceType,LMPHostType> &&
!std::is_same_v<DeviceType,LMPDeviceType>);
request->set_kokkos_device(std::is_same_v<DeviceType,LMPDeviceType>);
if (neighflag == FULL) request->enable_full();
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
template<class DeviceType>
double PairLJSPICACoulLongKokkos<DeviceType>::init_one(int i, int j)
{
double cutone = PairLJSPICACoulLong::init_one(i,j);
k_params.h_view(i,j).lj1 = lj1[i][j];
k_params.h_view(i,j).lj2 = lj2[i][j];
k_params.h_view(i,j).lj3 = lj3[i][j];
k_params.h_view(i,j).lj4 = lj4[i][j];
k_params.h_view(i,j).offset = offset[i][j];
k_params.h_view(i,j).cut_ljsq = cut_ljsq[i][j];
k_params.h_view(i,j).cut_coulsq = cut_coulsq;
k_params.h_view(i,j).lj_type = lj_type[i][j];
k_params.h_view(j,i) = k_params.h_view(i,j);
if (i<MAX_TYPES_STACKPARAMS+1 && j<MAX_TYPES_STACKPARAMS+1) {
m_params[i][j] = m_params[j][i] = k_params.h_view(i,j);
m_cutsq[j][i] = m_cutsq[i][j] = cutone*cutone;
m_cut_ljsq[j][i] = m_cut_ljsq[i][j] = cut_ljsq[i][j];
m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsq;
}
k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsq[i][j];
k_cutsq.template modify<LMPHostType>();
k_cut_ljsq.template modify<LMPHostType>();
k_params.template modify<LMPHostType>();
return cutone;
}
namespace LAMMPS_NS {
template class PairLJSPICACoulLongKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class PairLJSPICACoulLongKokkos<LMPHostType>;
#endif
}

148
src/KOKKOS/pair_lj_spica_coul_long_kokkos.h Normal file
View File

@ -0,0 +1,148 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(lj/spica/coul/long/kk,PairLJSPICACoulLongKokkos<LMPDeviceType>);
PairStyle(lj/spica/coul/long/kk/device,PairLJSPICACoulLongKokkos<LMPDeviceType>);
PairStyle(lj/spica/coul/long/kk/host,PairLJSPICACoulLongKokkos<LMPHostType>);
PairStyle(lj/sdk/coul/long/kk,PairLJSPICACoulLongKokkos<LMPDeviceType>);
PairStyle(lj/sdk/coul/long/kk/device,PairLJSPICACoulLongKokkos<LMPDeviceType>);
PairStyle(lj/sdk/coul/long/kk/host,PairLJSPICACoulLongKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_PAIR_LJ_SPICA_COUL_LONG_KOKKOS_H
#define LMP_PAIR_LJ_SPICA_COUL_LONG_KOKKOS_H
#include "pair_kokkos.h"
#include "pair_lj_spica_coul_long.h"
#include "neigh_list_kokkos.h"
namespace LAMMPS_NS {
template<class DeviceType>
class PairLJSPICACoulLongKokkos : public PairLJSPICACoulLong {
public:
enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
enum {COUL_FLAG=1};
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
PairLJSPICACoulLongKokkos(class LAMMPS *);
~PairLJSPICACoulLongKokkos() override;
void compute(int, int) override;
void init_tables(double cut_coul, double *cut_respa) override;
void init_style() override;
double init_one(int, int) override;
struct params_lj_spica_coul {
KOKKOS_INLINE_FUNCTION
params_lj_spica_coul() {cut_ljsq=0;cut_coulsq=0;lj1=0;lj2=0;lj3=0;lj4=0;offset=0;lj_type=0;};
KOKKOS_INLINE_FUNCTION
params_lj_spica_coul(int /*i*/) {cut_ljsq=0;cut_coulsq=0;lj1=0;lj2=0;lj3=0;lj4=0;offset=0;lj_type=0;};
F_FLOAT cut_ljsq,cut_coulsq,lj1,lj2,lj3,lj4,offset;
int lj_type;
};
protected:
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fpair(const F_FLOAT& rsq, const int& i, const int& j, const int& itype, const int& jtype) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_evdwl(const F_FLOAT& rsq, const int& i, const int& j, const int& itype, const int& jtype) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_fcoul(const F_FLOAT& rsq, const int& i, const int& j, const int& itype,
const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION
F_FLOAT compute_ecoul(const F_FLOAT& rsq, const int& i, const int& j, const int& itype,
const int& jtype, const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const;
Kokkos::DualView<params_lj_spica_coul**,Kokkos::LayoutRight,DeviceType> k_params;
typename Kokkos::DualView<params_lj_spica_coul**,Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
params_lj_spica_coul m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1]; // hardwired to space for 12 atom types
F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
typename AT::t_x_array_randomread x;
typename AT::t_x_array c_x;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread type;
typename AT::t_float_1d_randomread q;
DAT::tdual_efloat_1d k_eatom;
DAT::tdual_virial_array k_vatom;
typename AT::t_efloat_1d d_eatom;
typename AT::t_virial_array d_vatom;
int newton_pair;
typename AT::tdual_ffloat_2d k_cutsq, k_cut_ljsq;
typename AT::t_ffloat_2d d_cutsq, d_cut_ljsq, d_cut_coulsq;
typename AT::t_ffloat_1d_randomread
d_rtable, d_drtable, d_ftable, d_dftable,
d_ctable, d_dctable, d_etable, d_detable;
int neighflag;
int nlocal,nall,eflag,vflag;
double special_lj[4];
double special_coul[4];
double qqrd2e;
void allocate() override;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,true,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,true,1,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALF,true,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALFTHREAD,true,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,false,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,false,1,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALF,false,0,CoulLongTable<1>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALFTHREAD,false,0,CoulLongTable<1>>;
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,FULL,0,CoulLongTable<1>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,FULL,1,CoulLongTable<1>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,HALF,0,CoulLongTable<1>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,HALFTHREAD,0,CoulLongTable<1>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute<PairLJSPICACoulLongKokkos,CoulLongTable<1>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,true,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,true,1,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALF,true,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALFTHREAD,true,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,false,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,FULL,false,1,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALF,false,0,CoulLongTable<0>>;
friend struct PairComputeFunctor<PairLJSPICACoulLongKokkos,HALFTHREAD,false,0,CoulLongTable<0>>;
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,FULL,0,CoulLongTable<0>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,FULL,1,CoulLongTable<0>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,HALF,0,CoulLongTable<0>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute_neighlist<PairLJSPICACoulLongKokkos,HALFTHREAD,0,CoulLongTable<0>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend EV_FLOAT pair_compute<PairLJSPICACoulLongKokkos,CoulLongTable<0>>(PairLJSPICACoulLongKokkos*,NeighListKokkos<DeviceType>*);
friend void pair_virial_fdotr_compute<PairLJSPICACoulLongKokkos>(PairLJSPICACoulLongKokkos*);
};
}
#endif
#endif

18
src/KOKKOS/pair_lj_spica_kokkos.cpp
View File

@ -70,7 +70,6 @@ void PairLJSPICAKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
eflag = eflag_in; eflag = eflag_in;
vflag = vflag_in; vflag = vflag_in;
if (neighflag == FULL) no_virial_fdotr_compute = 1; if (neighflag == FULL) no_virial_fdotr_compute = 1;
ev_init(eflag,vflag,0); ev_init(eflag,vflag,0);
@ -108,6 +107,8 @@ void PairLJSPICAKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
// loop over neighbors of my atoms // loop over neighbors of my atoms
copymode = 1;
EV_FLOAT ev = pair_compute<PairLJSPICAKokkos<DeviceType>,void >(this,(NeighListKokkos<DeviceType>*)list); EV_FLOAT ev = pair_compute<PairLJSPICAKokkos<DeviceType>,void >(this,(NeighListKokkos<DeviceType>*)list);
if (eflag) eng_vdwl += ev.evdwl; if (eflag) eng_vdwl += ev.evdwl;
@ -132,8 +133,13 @@ void PairLJSPICAKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
if (vflag_fdotr) pair_virial_fdotr_compute(this); if (vflag_fdotr) pair_virial_fdotr_compute(this);
copymode = 0;
} }
/* ----------------------------------------------------------------------
compute pair force between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType> template<class DeviceType>
template<bool STACKPARAMS, class Specialisation> template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
@ -152,6 +158,10 @@ compute_fpair(const F_FLOAT &rsq, const int &, const int &, const int &itype, co
return a* ( lj_1*r6inv*b - lj_2 * r2inv); return a* ( lj_1*r6inv*b - lj_2 * r2inv);
} }
/* ----------------------------------------------------------------------
compute pair potential energy between atoms i and j
---------------------------------------------------------------------- */
template<class DeviceType> template<class DeviceType>
template<bool STACKPARAMS, class Specialisation> template<bool STACKPARAMS, class Specialisation>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
@ -166,18 +176,14 @@ compute_evdwl(const F_FLOAT &rsq, const int &, const int &, const int &itype, co
if (ljt == LJ12_4) { if (ljt == LJ12_4) {
const F_FLOAT r4inv=r2inv*r2inv; const F_FLOAT r4inv=r2inv*r2inv;
return r4inv*(lj_3*r4inv*r4inv - lj_4) - offset; return r4inv*(lj_3*r4inv*r4inv - lj_4) - offset;
} else if (ljt == LJ9_6) { } else if (ljt == LJ9_6) {
const F_FLOAT r3inv = r2inv*sqrt(r2inv); const F_FLOAT r3inv = r2inv*sqrt(r2inv);
const F_FLOAT r6inv = r3inv*r3inv; const F_FLOAT r6inv = r3inv*r3inv;
return r6inv*(lj_3*r3inv - lj_4) - offset; return r6inv*(lj_3*r3inv - lj_4) - offset;
} else if (ljt == LJ12_6) { } else if (ljt == LJ12_6) {
const double r6inv = r2inv*r2inv*r2inv; const double r6inv = r2inv*r2inv*r2inv;
return r6inv*(lj_3*r6inv - lj_4) - offset; return r6inv*(lj_3*r6inv - lj_4) - offset;
} else if (ljt == LJ12_5) { } else if (ljt == LJ12_5) {
const F_FLOAT r5inv = r2inv*r2inv*sqrt(r2inv); const F_FLOAT r5inv = r2inv*r2inv*sqrt(r2inv);
const F_FLOAT r7inv = r5inv*r2inv; const F_FLOAT r7inv = r5inv*r2inv;
@ -273,8 +279,6 @@ double PairLJSPICAKokkos<DeviceType>::init_one(int i, int j)
return cutone; return cutone;
} }
namespace LAMMPS_NS { namespace LAMMPS_NS {
template class PairLJSPICAKokkos<LMPDeviceType>; template class PairLJSPICAKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU #ifdef LMP_KOKKOS_GPU

9
unittest/force-styles/test_angle_style.cpp
View File

@ -645,8 +645,13 @@ TEST(AngleStyle, single)
"extra/angle/per/atom 2 extra/special/per/atom 2", "extra/angle/per/atom 2 extra/special/per/atom 2",
nangletypes)); nangletypes));
command("pair_style zero 8.0"); if (utils::strmatch(test_config.angle_style, "^spica")) {
command("pair_coeff * *"); command("pair_style lj/spica 8.0");
command("pair_coeff * * lj9_6 0.02 2.5");
} else {
command("pair_style zero 8.0");
command("pair_coeff * *");
}
command("angle_style " + test_config.angle_style); command("angle_style " + test_config.angle_style);
Angle *angle = lmp->force->angle; Angle *angle = lmp->force->angle;

91
unittest/force-styles/tests/angle-spica.yaml Normal file
View File

@ -0,0 +1,91 @@
---
lammps_version: 7 Feb 2024
tags:
date_generated: Wed Jun 5 18:31:11 2024
epsilon: 1e-12
skip_tests:
prerequisites: ! |
atom full
angle spica
pair lj/spica
pre_commands: ! |
variable write_data_pair index ij
post_commands: ! |
pair_style lj/spica 8.0
input_file: in.spica
angle_style: spica
angle_coeff: ! |
1 33.5 110.1
2 46.1 111.3
3 40.0 120.0
4 33.0 108.5
equilibrium: 4 1.9216075064457565 1.9425514574696887 2.0943951023931953 1.8936822384138474
extract: ! ""
natoms: 29
init_energy: 38.36438529349082
init_stress: ! |2-
6.2288484633192937e+01 -2.4958587357033732e+01 1.5261156310077459e+02 2.7527094178009044e+01 4.2708401447133369e+01 -2.5265950282815652e+01
init_forces: ! |2
1 1.7230431704651284e+01 4.0071311825737794e+01 2.5153895595391262e+01
2 1.4681450443715043e+01 9.4049099264163125e+00 -2.0102364558934152e+01
3 -1.9308548347980725e+01 -2.8460741684874360e+01 6.4339170989100403e+00
4 -5.2290857572683347e+00 -9.4927850504411957e+00 1.8003472602849664e+00
5 -1.6382493884699279e+01 -1.4784175400797803e+01 -6.1931262199840624e+00
6 -5.4493064751680564e+01 7.5613258287794565e+01 6.9941734306087866e+01
7 -4.7395486031142198e+01 3.7640428970968820e+00 -2.4566126298992822e+02
8 4.8960191044380039e+01 -4.5712855041323763e+01 1.9415466008820243e+02
9 -8.3924213860655783e+00 2.4590103192063562e+01 3.2878160380265854e+01
10 6.8038129570966305e+01 -4.6569300903873433e+01 -7.6671701712056688e+01
11 -2.1258622382772373e+01 -7.4008599257925383e+00 8.3610564111488799e+00
12 5.9901920466154710e+00 -3.0542621446562030e+00 -7.5713636171526382e+00
13 1.4189850824692691e+00 1.0835451791895949e+00 4.9993493214959894e-01
14 1.1355220599075533e+01 2.4997321466757656e+00 -5.2256561116849181e+00
15 1.3933035015335755e+00 -2.0786326583247088e+00 8.7283364302611144e+00
16 4.5875574605935235e+01 -2.7754880718666161e+01 -1.0582764215323260e+02
17 -4.2451226361036497e+01 2.8306085448648393e+01 1.1928362036932688e+02
18 5.2926497658556788e-02 7.0441277800577096e-01 -2.9584690206819020e+00
19 -8.6276212064612834e-01 -1.0756990560638664e+00 1.2826229074637665e+00
20 7.9148971075772434e-01 3.3735689913016514e-01 1.7102043999486722e+00
21 -5.9956840446092086e+00 -6.6440354071373493e+00 1.8013592125677949e+01
22 -1.4065918957775763e+01 -4.6031498182035513e+00 -1.5782854315475038e+01
23 2.0049558879496164e+01 1.1260925334020174e+01 -2.2177902350101499e+00
24 3.2508269160599186e+00 -1.9979416528015960e+01 1.0409971951734974e+01
25 -1.6060493526274090e+01 9.4218054960789366e-01 -1.2854339891002313e+01
26 1.2797730884335980e+01 1.9024610945429981e+01 2.4201678955447798e+00
27 5.3169168679178105e+00 -2.2776165356368182e+01 9.1726846500285717e+00
28 -1.9217628865768805e+01 7.5096140408582208e+00 -1.2798872429928149e+01
29 1.3910508062151619e+01 1.5274870243863951e+01 3.6205364526432242e+00
run_energy: 38.176417113761396
run_stress: ! |2-
6.1357382915550481e+01 -2.4715660959892247e+01 1.5186466302604836e+02 2.7218614858138846e+01 4.2112515489269541e+01 -2.5607593286714120e+01
run_forces: ! |2
1 1.7145760907927180e+01 4.0174672565571342e+01 2.5205552180237635e+01
2 1.4615170171901832e+01 9.2708448155131897e+00 -2.0078351161345097e+01
3 -1.9027637897644411e+01 -2.8642960432918144e+01 6.2314807077302952e+00
4 -5.2736618463889080e+00 -9.3424901289763778e+00 1.9422015848783882e+00
5 -1.6468832074580213e+01 -1.4729606520686309e+01 -6.2023333708615995e+00
6 -5.3718294227422540e+01 7.4801671388546964e+01 6.8368190594240389e+01
7 -4.7118098933035071e+01 3.8739803262288781e+00 -2.4290995583077535e+02
8 4.7876325968669121e+01 -4.4923173742504034e+01 1.9297993934786507e+02
9 -8.3824572561250399e+00 2.4542635562122975e+01 3.2839238074098198e+01
10 6.8198580866574645e+01 -4.6830380612887495e+01 -7.6629764788623447e+01
11 -2.1232237676690062e+01 -7.3683045541850962e+00 8.3643364257797792e+00
12 5.6908584692226185e+00 -2.9956405260312478e+00 -7.3282063658600993e+00
13 1.3684378256951288e+00 9.7105537426706034e-01 5.7801503997630133e-01
14 1.1379285236061559e+01 2.6704814459393837e+00 -5.2142720262734468e+00
15 1.5246096380310767e+00 -1.9886210527177606e+00 8.4529064022347313e+00
16 4.5766629638536287e+01 -2.7729192853260034e+01 -1.0555847936417963e+02
17 -4.2311907917867423e+01 2.8269504758036135e+01 1.1894203420097222e+02
18 4.5978604531469895e-02 6.0605282956542328e-01 -2.5489959663795330e+00
19 -7.4282409641952973e-01 -9.2821937232343465e-01 1.1099668747784261e+00
20 6.7850316938866917e-01 2.8824602749664086e-01 1.4733945401265087e+00
21 -6.1707359561298754e+00 -6.7761831890615607e+00 1.8428059527969836e+01
22 -1.3998960882672469e+01 -4.5334962854737153e+00 -1.5964259983712697e+01
23 2.0157648747423554e+01 1.1323428247127397e+01 -2.4508469770619676e+00
24 3.5283602170234043e+00 -2.0588385440875950e+01 1.0846305899711011e+01
25 -1.6431076320543696e+01 1.1358029038847810e+00 -1.3239288833419160e+01
26 1.2890779580326928e+01 1.9439960296533382e+01 2.3687830191265684e+00
27 5.4188347979018312e+00 -2.2961492507008920e+01 9.2221110477962096e+00
28 -1.9337322641186159e+01 7.5984768842531718e+00 -1.2868573316667142e+01
29 1.3928283887490117e+01 1.5371333793823325e+01 3.6408125176375652e+00
...

233
unittest/force-styles/tests/data.spica Normal file
View File

@ -0,0 +1,233 @@
LAMMPS data file via write_data, version 5 May 2020, timestep = 0
29 atoms
5 atom types
24 bonds
5 bond types
30 angles
4 angle types
31 dihedrals
5 dihedral types
2 impropers
2 improper types
-6.024572 8.975428 xlo xhi
-7.692866 7.307134 ylo yhi
-8.086924 6.913076 zlo zhi
Masses
1 12.0107
2 4.00794
3 14.0067
4 15.9994
5 15.9994
PairIJ Coeffs # lj/spica
1 1 lj9_6 0.02 2.5
1 2 lj9_6 0.01 1.58114
1 3 lj9_6 0.02 2.82843
1 4 lj9_6 0.0173205 2.78388
1 5 lj9_6 0.0173205 2.78388
2 2 lj12_4 0.005 1.0
2 3 lj12_4 0.01 1.78885
2 4 lj12_4 0.005 0.5
2 5 lj12_4 0.00866025 1.76068 8
3 3 lj12_6 0.02 3.2 8
3 4 lj12_6 0.0173205 3.1496 8
3 5 lj12_6 0.0173205 3.1496 8
4 4 lj9_6 0.015 3.1 8
4 5 lj9_6 0.015 3.1 8
5 5 lj9_6 0.015 3.1 8
Bond Coeffs # zero
1 1.5
2 1.1
3 1.3
4 1.2
5 1
Angle Coeffs # spica
1 33.5 110.1
2 46.1 111.3
3 40 120
4 33 108.5
Dihedral Coeffs # zero
1
2
3
4
5
Improper Coeffs # zero
1
2
Atoms # full
10 2 1 7.0000000000000007e-02 2.0185283555536988e+00 -1.4283966846517357e+00 -9.6733527271133024e-01 0 0 0
11 2 2 8.9999999999999997e-02 1.7929780509347666e+00 -1.9871047540768743e+00 -1.8840626643185674e+00 0 0 0
12 2 1 -2.7000000000000002e-01 3.0030247876861225e+00 -4.8923319967572748e-01 -1.6188658531537248e+00 0 0 0
13 2 2 8.9999999999999997e-02 4.0447273787895934e+00 -9.0131998547446246e-01 -1.6384447268320836e+00 0 0 0
14 2 2 8.9999999999999997e-02 2.6033152817257075e+00 -4.0789761505963579e-01 -2.6554413538823063e+00 0 0 0
2 1 2 3.1000000000000000e-01 3.0197083955402204e-01 2.9515239068888608e+00 -8.5689735572907566e-01 0 0 0
3 1 1 -2.0000000000000000e-02 -6.9435377880558602e-01 1.2440473127136711e+00 -6.2233801468892025e-01 0 0 0
4 1 2 8.9999999999999997e-02 -1.5771614164685133e+00 1.4915333140468066e+00 -1.2487126845040522e+00 0 0 0
6 1 1 5.1000000000000001e-01 2.9412607937706009e-01 2.2719282656652909e-01 -1.2843094067857870e+00 0 0 0
7 1 4 -5.1000000000000001e-01 3.4019871062879609e-01 -9.1277350075786561e-03 -2.4633113224304561e+00 0 0 0
19 3 2 4.2359999999999998e-01 1.5349125211132961e+00 2.6315969880333707e+00 -4.2472859440220647e+00 0 0 0
15 2 2 8.9999999999999997e-02 2.9756315249791303e+00 5.6334269722969288e-01 -1.2437650754599008e+00 0 0 0
18 3 4 -8.4719999999999995e-01 2.1384791188033843e+00 3.0177261773770208e+00 -3.5160827596876225e+00 0 0 0
20 3 2 4.2359999999999998e-01 2.7641167828863153e+00 3.6833419064000221e+00 -3.9380850623312638e+00 0 0 0
8 2 3 -4.6999999999999997e-01 1.1641187171852805e+00 -4.8375305955385234e-01 -6.7659823767368688e-01 0 0 0
9 2 2 3.1000000000000000e-01 1.3777459838125838e+00 -2.5366338669522998e-01 2.6877644730326306e-01 0 0 0
16 2 1 5.1000000000000001e-01 2.6517554244980306e+00 -2.3957110424978438e+00 3.2908335999178327e-02 0 0 0
17 2 4 -5.1000000000000001e-01 2.2309964792710639e+00 -2.1022918943319384e+00 1.1491948328949437e+00 0 0 0
1 1 3 -4.6999999999999997e-01 -2.7993683669226832e-01 2.4726588069312840e+00 -1.7200860244148433e-01 0 0 0
5 1 2 8.9999999999999997e-02 -8.9501761359359255e-01 9.3568128743071344e-01 4.0227731871484346e-01 0 0 0
21 4 5 -8.4719999999999995e-01 4.9064454390208301e+00 -4.0751205255383196e+00 -3.6215576073601046e+00 0 0 0
22 4 2 4.2359999999999998e-01 4.3687453488627543e+00 -4.2054270536772504e+00 -4.4651491269372565e+00 0 0 0
23 4 2 4.2359999999999998e-01 5.7374928154769504e+00 -3.5763355905184966e+00 -3.8820297194230728e+00 0 0 0
24 5 5 -8.4719999999999995e-01 2.0684115301174013e+00 3.1518221747664397e+00 3.1554242678474576e+00 0 0 0
25 5 2 4.2359999999999998e-01 1.2998381073113014e+00 3.2755513587518097e+00 2.5092990173114837e+00 0 0 0
26 5 2 4.2359999999999998e-01 2.5807438597688113e+00 4.0120175892854135e+00 3.2133398379059099e+00 0 0 0
27 6 5 -8.4719999999999995e-01 -1.9613581876744359e+00 -4.3556300596085160e+00 2.1101467673534788e+00 0 0 0
28 6 2 4.2359999999999998e-01 -2.7406520384725965e+00 -4.0207251278130975e+00 1.5828689861678511e+00 0 0 0
29 6 2 4.2359999999999998e-01 -1.3108232656499081e+00 -3.5992986322410760e+00 2.2680459788743503e+00 0 0 0
Velocities
1 7.7867804888392077e-04 5.8970331623292821e-04 -2.2179517633030531e-04
2 2.7129529964126462e-03 4.6286427111164284e-03 3.5805549693846352e-03
3 -1.2736791029204805e-03 1.6108674226414498e-03 -3.3618185901550799e-04
4 -9.2828595122009308e-04 -1.2537885319521818e-03 -4.1204974953432108e-03
5 -1.1800848061603740e-03 7.5424401975844038e-04 6.9023177964912290e-05
6 -3.0914004879905335e-04 1.2755385764678133e-03 7.9574303350202582e-04
7 -1.1037894966874103e-04 -7.6764845099077425e-04 -7.7217630460203659e-04
8 3.9060281273221989e-04 -8.1444231918053418e-04 1.5134641148324972e-04
9 1.2475530960659720e-03 -2.6608454451432528e-03 1.1117602907112732e-03
10 4.5008983776042893e-04 4.9530197647538077e-04 -2.3336234361093645e-04
11 -3.6977669078869707e-04 -1.5289071951960539e-03 -2.9176389881837113e-03
12 1.0850834530183159e-03 -6.4965897903201833e-04 -1.2971152622619948e-03
13 4.0754559196230639e-03 3.5043502394946119e-03 -7.8324487687854666e-04
14 -1.3837220448746613e-04 -4.0656048637594394e-03 -3.9333461173944500e-03
15 -4.3301707382721859e-03 -3.1802661664634938e-03 3.2037919043360571e-03
16 -9.6715751018414326e-05 -5.0016572678960377e-04 1.4945658875149626e-03
17 6.5692180538157174e-04 3.6635779995305095e-04 8.3495414466050911e-04
18 -6.0936815808025862e-04 -9.3774557532468582e-04 -3.3558072507805731e-04
19 -6.9919768291957119e-04 -3.6060777270430031e-03 4.2833405289822791e-03
20 4.7777805013736515e-03 5.1003745845520452e-03 1.8002873923729241e-03
21 -9.5568188553430398e-04 1.6594630943762931e-04 -1.8199788009966615e-04
22 -3.3137518957653462e-03 -2.8683968287936054e-03 3.6384389958326871e-03
23 2.4209481134686401e-04 -4.5457709985051130e-03 2.7663581642115042e-03
24 2.5447450568861086e-04 4.8412447786110117e-04 -4.8021914527341357e-04
25 4.3722771097312743e-03 -4.5184411669545515e-03 2.5200952006556795e-03
26 -1.9250110555001179e-03 -3.0342169883610837e-03 3.5062814567984532e-03
27 -2.6510179146429716e-04 3.6306203629019116e-04 -5.6235585400647747e-04
28 -2.3068708109787484e-04 -8.5663070212203200e-04 2.1302563179109169e-03
29 -2.5054744388303732e-03 -1.6773997805290820e-04 2.8436699761004796e-03
Bonds
1 5 1 2
2 3 1 3
3 2 3 4
4 2 3 5
5 1 3 6
6 3 6 8
7 4 6 7
8 5 8 9
9 3 8 10
10 2 10 11
11 1 10 12
12 1 10 16
13 2 12 13
14 2 12 14
15 2 12 15
16 4 16 17
17 5 18 19
18 5 18 20
19 5 21 22
20 5 21 23
21 5 24 25
22 5 24 26
23 5 27 28
24 5 27 29
Angles
1 4 2 1 3
2 4 1 3 5
3 4 1 3 4
4 4 1 3 6
5 4 4 3 5
6 2 5 3 6
7 2 4 3 6
8 3 3 6 7
9 3 3 6 8
10 3 7 6 8
11 2 6 8 9
12 2 9 8 10
13 3 6 8 10
14 2 8 10 11
15 3 8 10 16
16 2 11 10 12
17 1 12 10 16
18 1 8 10 12
19 2 11 10 16
20 2 10 12 15
21 2 10 12 14
22 2 10 12 13
23 4 13 12 15
24 4 13 12 14
25 4 14 12 15
26 4 10 16 17
27 1 19 18 20
28 1 22 21 23
29 1 25 24 26
30 1 28 27 29
Dihedrals
1 2 2 1 3 6
2 2 2 1 3 4
3 3 2 1 3 5
4 1 1 3 6 8
5 1 1 3 6 7
6 5 4 3 6 8
7 5 4 3 6 7
8 5 5 3 6 8
9 5 5 3 6 7
10 4 3 6 8 9
11 3 3 6 8 10
12 3 7 6 8 9
13 4 7 6 8 10
14 2 6 8 10 12
15 2 6 8 10 16
16 2 6 8 10 11
17 2 9 8 10 12
18 4 9 8 10 16
19 5 9 8 10 11
20 5 8 10 12 13
21 1 8 10 12 14
22 5 8 10 12 15
23 4 8 10 16 17
24 5 11 10 12 13
25 5 11 10 12 14
26 5 11 10 12 15
27 2 11 10 16 17
28 2 12 10 16 17
29 5 16 10 12 13
30 5 16 10 12 14
31 5 16 10 12 15
Impropers
1 1 6 3 8 7
2 2 8 6 10 9

30
unittest/force-styles/tests/in.spica Normal file
View File

@ -0,0 +1,30 @@
variable newton_pair index on
variable newton_bond index on
variable bond_factor index 0.10
variable angle_factor index 0.25
variable dihedral_factor index 0.50
variable units index real
variable input_dir index .
variable data_file index ${input_dir}/data.spica
variable pair_style index 'lj/spica 8.0'
variable bond_style index zero
variable angle_style index spica
variable dihedral_style index zero
variable improper_style index zero
variable t_target index 100.0
atom_style full
atom_modify map array
neigh_modify delay 2 every 2 check no
units ${units}
timestep 0.1
newton ${newton_pair} ${newton_bond}
special_bonds lj/coul ${bond_factor} ${angle_factor} ${dihedral_factor}
pair_style ${pair_style}
bond_style ${bond_style}
angle_style ${angle_style}
dihedral_style ${dihedral_style}
improper_style ${improper_style}
read_data ${data_file}