diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index aaa899ef66..c2753b2afc 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -104,6 +104,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`pe/tally ` * :doc:`plasticity/atom ` * :doc:`pressure ` + * :doc:`pressure/alchemy ` * :doc:`pressure/uef ` * :doc:`property/atom ` * :doc:`property/chunk ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index d275b33eba..5b331af1aa 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -29,6 +29,7 @@ OPT. * :doc:`adapt/fep ` * :doc:`addforce ` * :doc:`addtorque ` + * :doc:`alchemy ` * :doc:`amoeba/bitorsion ` * :doc:`amoeba/pitorsion ` * :doc:`append/atoms ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 72ea31bbb4..3d74b3884e 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -258,6 +258,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`pe/tally ` - potential energy between two groups of atoms via the tally callback mechanism * :doc:`plasticity/atom ` - Peridynamic plasticity for each atom * :doc:`pressure ` - total pressure and pressure tensor +* :doc:`pressure/alchemy ` - mixed system total pressure and pressure tensor for :doc:`fix alchemy ` runs * :doc:`pressure/uef ` - pressure tensor in the reference frame of an applied flow field * :doc:`property/atom ` - convert atom attributes to per-atom vectors/arrays * :doc:`property/chunk ` - extract various per-chunk attributes diff --git a/doc/src/compute_pressure_alchemy.rst b/doc/src/compute_pressure_alchemy.rst new file mode 100644 index 0000000000..bdf9802e20 --- /dev/null +++ b/doc/src/compute_pressure_alchemy.rst @@ -0,0 +1,80 @@ +.. index:: compute pressure/alchemy + +compute pressure/alchemy command +================================ + +Syntax +"""""" + +.. code-block:: LAMMPS + + compute ID group-ID pressure/alchemy fix-ID + +* ID, group-ID are documented in :doc:`compute ` command +* pressure/alchemy = style name of this compute command +* fix-ID = ID of :doc:`fix alchemy ` command + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix trans all alchemy + compute mixed all pressure/alchemy trans + thermo_modify press mixed + +Description +""""""""""" + +.. versionadded:: TBD + +Define a compute style that makes the "mixed" system pressure available +for a system that uses the :doc:`fix alchemy ` command to +transform one topology to another. This can be used in combination with +either :doc:`thermo_modify press ` or :doc:`fix_modify +press ` to output and access a pressure consistent with the +simulated combined two topology system. + +The actual pressure is determined with :doc:`compute pressure +` commands that are internally used by :doc:`fix +alchemy ` for each topology individually and then combined. +This command just extracts the information from the fix. + +The ``examples/PACKAGES/alchemy`` folder contains an example input for this command. + +---------- + +Output info +""""""""""" + +This compute calculates a global scalar (the pressure) and a global +vector of length 6 (the pressure tensor), which can be accessed by +indices 1--6. These values can be used by any command that uses global +scalar or vector values from a compute as input. See the :doc:`Howto +output ` page for an overview of LAMMPS output options. + +The ordering of values in the symmetric pressure tensor is as follows: +:math:`p_{xx},` :math:`p_{yy},` :math:`p_{zz},` :math:`p_{xy},` +:math:`p_{xz},` :math:`p_{yz}.` + +The scalar and vector values calculated by this compute are "intensive". +The scalar and vector values will be in pressure :doc:`units `. + +Restrictions +"""""""""""" + +This compute is part of the REPLICA package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + + +Related commands +"""""""""""""""" + +:doc:`fix alchemy `, :doc:`compute pressure `, +:doc:`thermo_modify `, :doc:`fix_modify ` + +Default +""""""" + +none diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 2dfe97a3ec..71f8bf9e16 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -181,6 +181,7 @@ accelerated styles exist. * :doc:`adapt/fep ` - enhanced version of fix adapt * :doc:`addforce ` - add a force to each atom * :doc:`addtorque ` - add a torque to a group of atoms +* :doc:`alchemy ` - perform an "alchemical transformation" between two partitions * :doc:`amoeba/bitorsion ` - torsion/torsion terms in AMOEBA force field * :doc:`amoeba/pitorsion ` - 6-body terms in AMOEBA force field * :doc:`append/atoms ` - append atoms to a running simulation diff --git a/doc/src/fix_alchemy.rst b/doc/src/fix_alchemy.rst new file mode 100644 index 0000000000..90c436139e --- /dev/null +++ b/doc/src/fix_alchemy.rst @@ -0,0 +1,137 @@ +.. index:: fix alchemy + +fix alchemy command +=================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID alchemy v_name + +* ID, group-ID are documented in :doc:`fix ` command +* alchemy = style name of this fix command +* v_name = variable with name that determines the :math:`\lambda_p` value + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix trans all alchemy v_ramp + +Description +""""""""""" + +.. versionadded:: TBD + +This fix command enables running an "alchemical transformation" MD +simulation between two topologies (i.e. the same number and positions of +atoms, but differences in atom parameters like type, charge, bonds, +angles and so on). For this a :ref:`multi-partition run ` is +required with exactly two partitions. During the MD run, the fix will +will determine a factor, :math:`\lambda_p`, for each partition *p* that +will be taken from an equal style or equivalent :doc:`variable +`. Typically, this variable would be chose to linearly ramp +*down* from 1.0 to 0.0 for the *first* partition (*p=0*) and linearly +ramp *up* from 0.0 to 1.0 for the *second* partition (*p=1*). The +forces used for the propagation of the atoms will be the sum of the +forces of the two systems combined and scaled with their respective +:math:`\lambda_p` factor. This allows to perform transformations that +are not easily possible with :doc:`pair style hybrid/scaled +`, :doc:`fix adapt ` or :doc:`fix adapt/fep +`. + +.. note:: + + Since the definition of the variable to provide the :math:`\lambda_p` is + independent in the two partitions, no check is made that the two values + remain between 0.0 and 1.0 and that they add up to 1.0. So care needs to + be taken when defining those variables that this is the case. + +Due to the specifics of the implementation, the initial geometry and +dimensions of the system must be exactly the same and the fix will +synchronize them during the run. It is thus not possible to initialize +the two partitions by reading different data files or creating different +systems from scratch, but rather they have to be started from the same +system and then the desired modifications need to be applied to the +system of the second partition. The commands :doc:`pair style +hybrid/scaled `, :doc:`fix adapt ` or :doc:`fix +adapt/fep ` could be used for simulations where the +requirements for fix alchemy are not given. + +The commands below demonstrate how the setup for the second partition +can be done for the example of transforming a pure copper system into a +copper/aluminum bronze. + +.. code-block:: LAMMPS + + variable name world pure alloy + + create_box 2 box + create_atoms 1 box + pair_style eam/alloy + pair_coeff * * AlCu.eam.alloy Cu Al + + # replace 5% of copper with aluminum on the second partition only + variable name world pure alloy + if "${name} == alloy" then & + "set type 1 type/fraction 2 0.05 6745234" + + # define ramp variable to combine the two different partitions + if "${name} == pure" then & + "variable ramp equal ramp(1.0,0.0)" & + else & + "variable ramp equal ramp(0.0,1.0)" + + fix 2 all alchemy v_ramp + + +The ``examples/PACKAGES/alchemy`` folder contains complete example +inputs for this command. + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files `. +None of the :doc:`fix_modify ` options are relevant to this fix. + +This fix stores a global scalar (the current value of :math:`\lambda_p`) +and a global vector of length 3 which contains the potential energy of +the first partition, the second partition and the combined value, +respectively. The global scalar is unitless and "intensive", the vector +is in :doc:`energy units ` and "extensive". These values can be +used by any command that uses a global value from a fix as input. See +the :doc:`Howto output ` doc page for an overview of +LAMMPS output options. + +This fix is not invoked during :doc:`energy minimization `. + +Restrictions +"""""""""""" + +This fix is part of the REPLICA package. It is only enabled if LAMMPS +was built with that package. See the :doc:`Build package +` page for more info. + +There may be only one instance of this fix in use at any time. + +This fix requires to perform a :ref:`multi-partition run ` +with *exactly* two partitions. + +This fix is *not* compatible with :doc:`load balancing `. + +Related commands +"""""""""""""""" + +:doc:`compute pressure/alchemy ` command, +:doc:`fix adapt ` command, :doc:`fix adapt/fep ` +command, :doc:`pair_style hybrid/scaled ` command. + +Default +""""""" + +none diff --git a/examples/PACKAGES/alchemy/AlCu.eam.alloy b/examples/PACKAGES/alchemy/AlCu.eam.alloy new file mode 120000 index 0000000000..c9470b2084 --- /dev/null +++ b/examples/PACKAGES/alchemy/AlCu.eam.alloy @@ -0,0 +1 @@ +../../../potentials/AlCu.eam.alloy \ No newline at end of file diff --git a/examples/PACKAGES/alchemy/h2o.mol b/examples/PACKAGES/alchemy/h2o.mol new file mode 100644 index 0000000000..e5a5e4fe93 --- /dev/null +++ b/examples/PACKAGES/alchemy/h2o.mol @@ -0,0 +1,62 @@ +# Water molecule. SPC/E model. + +3 atoms +2 bonds +1 angles + +Coords + +1 1.12456 0.09298 1.27452 +2 1.53683 0.75606 1.89928 +3 0.49482 0.56390 0.65678 + +Types + +1 1 +2 2 +3 2 + +Charges + +1 -0.8472 +2 0.4236 +3 0.4236 + +Bonds + +1 1 1 2 +2 1 1 3 + +Angles + +1 1 2 1 3 + +Shake Flags + +1 1 +2 1 +3 1 + +Shake Atoms + +1 1 2 3 +2 1 2 3 +3 1 2 3 + +Shake Bond Types + +1 1 1 1 +2 1 1 1 +3 1 1 1 + +Special Bond Counts + +1 2 0 0 +2 1 1 0 +3 1 1 0 + +Special Bonds + +1 2 3 +2 1 3 +3 1 2 diff --git a/examples/PACKAGES/alchemy/in.alloy b/examples/PACKAGES/alchemy/in.alloy new file mode 100644 index 0000000000..8cadcd7f39 --- /dev/null +++ b/examples/PACKAGES/alchemy/in.alloy @@ -0,0 +1,45 @@ +# Example for an alchemical transformation simulation of pure copper into an aluminium bronze + +# set up different names for two partitions +variable name world pure alloy + +# create pure copper system +units metal +lattice fcc 3.75 +region box block 0 6 0 6 0 6 +create_box 2 box + +timestep 0.002 +create_atoms 1 box +displace_atoms all random 0.3 0.3 0.3 57845645 +pair_style eam/alloy +pair_coeff * * AlCu.eam.alloy Cu Al + +# replace 5% of copper with aluminium on the second partition only +if "${name} == alloy" then & + "set type 1 type/fraction 2 0.05 6745234" + +# define ramp variable to combine the two different partitions +if "${name} == pure" then & + "variable ramp equal ramp(1.0,0.0)" & +else & + "variable ramp equal ramp(0.0,1.0)" + +velocity all create 5000.0 6567345 +fix 1 all nvt temp 2500.0 500.0 0.002 +fix 2 all alchemy v_ramp + +compute pressure all pressure/alchemy 2 + +# only need to output a dump file from one partition +#if "${name} == alloy" then & +# "dump 1 all custom 100 ${name}.lammpstrj id type element x y z" & +# "dump_modify 1 element Cu Al" + +thermo_style custom step temp press etotal density pe ke f_2 f_2[3] +thermo_modify colname f_2 lambda colname f_2[3] EPot_mixed +thermo_modify press pressure + +thermo 100 +run 10000 post no + diff --git a/examples/PACKAGES/alchemy/in.twowater b/examples/PACKAGES/alchemy/in.twowater new file mode 100644 index 0000000000..d39b0d9aaa --- /dev/null +++ b/examples/PACKAGES/alchemy/in.twowater @@ -0,0 +1,100 @@ +# Example for an alchemical transformation of two water molecules into a hydronium and hydroxyl ion +# WARNING: This input is intended for demonstrating the method only, +# the force field parameters are mostly made up and NOT suitable for production simulations. + +# set up different names for two partitions +variable name world twowater twoions + +units real +atom_style full +atom_modify map array +region box block -5 5 -5 5 -5 5 +boundary p p p +create_box 2 box bond/types 2 angle/types 2 & + extra/bond/per/atom 3 extra/angle/per/atom 3 extra/special/per/atom 3 + +mass 1 15.9994 +mass 2 1.008 + +pair_style lj/cut/coul/cut 10.0 +pair_coeff 1 1 0.1553 3.166 +pair_coeff 1 2 0.0 1.0 +pair_coeff 2 2 0.0 1.0 + +bond_style harmonic +bond_coeff * 1000.0 1.0 + +angle_style harmonic +angle_coeff * 100.0 109.47 + +molecule water h2o.mol + +# create the two molecules we want to transform ... +create_atoms 0 single -2.0 0.0 0.0 mol water 453624 +create_atoms 0 single 2.0 0.0 0.0 mol water 767353 + +# ... and put them in a group +group transform id 1:6 + +# now fill the rest of the box with more water +create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33 + +# change topology and settings for the two states +# we cannot simply create a different topology directly or +# load a different data file because the order and position +# of all atoms must be maintained across both replica + +# we first have to remove all topology data in the transform group +delete_bonds transform bond 1 +delete_bonds transform angle 1 remove + +# then generate different topologies for the two partitions. select by name. +if "${name} == twowater" then & + "create_bonds single/bond 2 1 2" & + "create_bonds single/bond 2 1 3" & + "create_bonds single/bond 2 4 5" & + "create_bonds single/bond 2 4 6" & + "create_bonds single/angle 2 2 1 3" & + "create_bonds single/angle 2 5 4 6" & +else & + "create_bonds single/bond 2 1 2" & + "create_bonds single/bond 2 3 4" & + "create_bonds single/bond 2 4 5" & + "create_bonds single/bond 2 4 6" & + "create_bonds single/angle 2 3 4 5" & + "create_bonds single/angle 2 5 4 6" & + "create_bonds single/angle 2 3 4 6" & + "set atom 1 charge -1.1354" & + "set atom 2 charge 0.1354" & + "set atom 3 charge 0.56775" & + "set atom 4 charge -0.70305" & + "set atom 5*6 charge 0.56775" + +velocity all create 300.0 5463576 +timestep 0.2 + +# define ramp variable to combine the two different partitions +if "${name} == twowater" then & + "variable ramp equal ramp(1.0,0.0)" & +else & + "variable ramp equal ramp(0.0,1.0)" + +# since the trajectory and forces are kept identical through fix alchemy, +# we can do fix npt simulations, but we must use the "mixed" pressure + +fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0 +fix transform all alchemy v_ramp +compute pressure all pressure/alchemy transform +fix_modify integrate press pressure + +# only need to output a dump file from one partition +# if "${name} == twowater" then & +# "dump 1 all atom 100 ${name}.lammpstrj" & +# "dump_modify 1 sort id" + +thermo_style custom step temp press etotal density pe ke f_transform f_transform[3] +thermo_modify colname f_transform lambda colname f_transform[3] EPot_mixed +thermo_modify press pressure + +thermo 100 +run 20000 diff --git a/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4 b/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4 new file mode 100644 index 0000000000..18375012d0 --- /dev/null +++ b/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4 @@ -0,0 +1,13 @@ +LAMMPS (8 Feb 2023) +Running on 2 partitions of processors +Starting alchemical transformation at 0% + Alchemical transformation progress: 10% + Alchemical transformation progress: 20% + Alchemical transformation progress: 30% + Alchemical transformation progress: 40% + Alchemical transformation progress: 50% + Alchemical transformation progress: 60% + Alchemical transformation progress: 70% + Alchemical transformation progress: 80% + Alchemical transformation progress: 90% + Alchemical transformation progress: 100% diff --git a/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.0 b/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.0 new file mode 100644 index 0000000000..c5538c3c3d --- /dev/null +++ b/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.0 @@ -0,0 +1,169 @@ +LAMMPS (8 Feb 2023) +Processor partition = 0 + using 1 OpenMP thread(s) per MPI task +# Example for an alchemical transformation simulation of pure copper into an aluminium bronze + +# set up different names for two partitions +variable name world pure alloy + +# create pure copper system +units metal +lattice fcc 3.75 +Lattice spacing in x,y,z = 3.75 3.75 3.75 +region box block 0 6 0 6 0 6 +create_box 2 box +Created orthogonal box = (0 0 0) to (22.5 22.5 22.5) + 1 by 1 by 2 MPI processor grid + +timestep 0.002 +create_atoms 1 box +Created 864 atoms + using lattice units in orthogonal box = (0 0 0) to (22.5 22.5 22.5) + create_atoms CPU = 0.001 seconds +displace_atoms all random 0.3 0.3 0.3 57845645 +Displacing atoms ... +pair_style eam/alloy +pair_coeff * * AlCu.eam.alloy Cu Al +Reading eam/alloy potential file AlCu.eam.alloy with DATE: 2008-10-01 + +# replace 5% of copper with aluminium on the second partition only +if "${name} == alloy" then "set type 1 type/fraction 2 0.05 6745234" + +# define ramp variable to combine the two different partitions +if "${name} == pure" then "variable ramp equal ramp(1.0,0.0)" else "variable ramp equal ramp(0.0,1.0)" +variable ramp equal ramp(1.0,0.0) + +velocity all create 5000.0 6567345 +fix 1 all nvt temp 2500.0 500.0 0.002 +fix 2 all alchemy v_ramp + +compute pressure all pressure/alchemy 2 + +# only need to output a dump file from one partition +#if "${name} == alloy" then # "dump 1 all custom 100 ${name}.lammpstrj id type element x y z" # "dump_modify 1 element Cu Al" + +thermo_style custom step temp press etotal density pe ke f_2 f_2[3] +thermo_modify colname f_2 lambda colname f_2[3] EPot_mixed +thermo_modify press pressure + +thermo 100 +run 10000 post no +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 8.6825 + ghost atom cutoff = 8.6825 + binsize = 4.34125, bins = 6 6 6 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair eam/alloy, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 3.474 | 3.474 | 3.474 Mbytes + Step Temp Press TotEng Density PotEng KinEng lambda EPot_mixed + 0 5000 2135632.8 10911.637 8.0039334 10353.88 557.75753 1 10353.88 + 100 2424.8703 98461.445 -2402.724 8.0039334 -2673.222 270.49793 0.99 -2673.3837 + 200 2404.4922 111478.55 -2381.5114 8.0039334 -2649.7361 268.22473 0.98 -2650.057 + 300 2388.2899 109429.75 -2387.7985 8.0039334 -2654.2158 266.41733 0.97 -2654.7075 + 400 2377.6418 97114.07 -2409.7114 8.0039334 -2674.9409 265.22953 0.96 -2675.755 + 500 2371.5028 98285.746 -2408.3248 8.0039334 -2672.8695 264.54471 0.95 -2673.6895 + 600 2368.0363 101125.13 -2402.9816 8.0039334 -2667.1396 264.15802 0.94 -2668.238 + 700 2364.8386 96916.491 -2410.1795 8.0039334 -2673.9808 263.80131 0.93 -2675.3394 + 800 2361.5815 101050.97 -2404.7579 8.0039334 -2668.1958 263.43797 0.92 -2669.6797 + 900 2355.1085 95290.542 -2416.4508 8.0039334 -2679.1667 262.7159 0.91 -2681.0381 + 1000 2344.579 90348.199 -2430.3553 8.0039334 -2691.8966 261.54132 0.9 -2693.7969 + 1100 2328.7808 87089.44 -2435.1426 8.0039334 -2694.9216 259.779 0.89 -2696.9199 + 1200 2310.0924 89628.642 -2438.4943 8.0039334 -2696.1886 257.69429 0.88 -2698.4983 + 1300 2288.2762 92771.824 -2430.0042 8.0039334 -2685.2648 255.26065 0.87 -2687.6406 + 1400 2258.3028 83029.221 -2451.9432 8.0039334 -2703.8603 251.91707 0.86 -2706.638 + 1500 2226.5131 83386 -2452.9658 8.0039334 -2701.3367 248.37089 0.85 -2704.0915 + 1600 2192.5777 81801.467 -2460.6692 8.0039334 -2705.2546 244.58535 0.84 -2708.3304 + 1700 2158.7421 84933.294 -2456.2782 8.0039334 -2697.0891 240.81093 0.83 -2700.142 + 1800 2123.8912 76700.779 -2474.2534 8.0039334 -2711.1767 236.92327 0.82 -2714.6978 + 1900 2091.3912 82568.985 -2469.6316 8.0039334 -2702.9295 233.29783 0.81 -2706.5643 + 2000 2061.868 79024.983 -2476.3262 8.0039334 -2706.3307 230.00448 0.8 -2710.7 + 2100 2036.179 72408.502 -2493.6936 8.0039334 -2720.8324 227.13884 0.79 -2725.1617 + 2200 2014.0161 81172.694 -2479.0525 8.0039334 -2703.719 224.66652 0.78 -2708.4618 + 2300 1996.2707 78029.129 -2485.3835 8.0039334 -2708.0706 222.687 0.77 -2712.5894 + 2400 1982.6409 79691.571 -2483.2441 8.0039334 -2704.4107 221.16658 0.76 -2709.4434 + 2500 1971.6774 75015.359 -2493.8754 8.0039334 -2713.819 219.94358 0.75 -2719.4159 + 2600 1963.0369 74845.99 -2493.3203 8.0039334 -2712.3 218.97972 0.74 -2718.1619 + 2700 1955.2025 76756.042 -2491.6282 8.0039334 -2709.7339 218.10578 0.73 -2716.0963 + 2800 1947.4326 70627.84 -2507.5853 8.0039334 -2724.8243 217.23904 0.72 -2730.7999 + 2900 1939.4308 70852.114 -2508.8277 8.0039334 -2725.1742 216.34642 0.71 -2731.5024 + 3000 1929.3184 75867.616 -2497.0826 8.0039334 -2712.301 215.21837 0.7 -2718.4146 + 3100 1916.202 68984.855 -2517.8316 8.0039334 -2731.5868 213.75522 0.69 -2738.0966 + 3200 1901.363 73304.621 -2508.0281 8.0039334 -2720.128 212.09991 0.68 -2727.1933 + 3300 1881.196 69447.29 -2513.9782 8.0039334 -2723.8285 209.85024 0.67 -2730.768 + 3400 1856.1833 64973.868 -2526.8095 8.0039334 -2733.8695 207.06004 0.66 -2741.4173 + 3500 1829.7246 64688.566 -2528.8664 8.0039334 -2732.9749 204.10853 0.65 -2740.986 + 3600 1800.1245 56451.766 -2547.6377 8.0039334 -2748.4443 200.80659 0.64 -2756.4626 + 3700 1771.5886 68268.147 -2530.5906 8.0039334 -2728.214 197.62337 0.63 -2735.8618 + 3800 1740.3894 64187.24 -2539.0772 8.0039334 -2733.2202 194.14306 0.62 -2742.67 + 3900 1709.6971 61697.296 -2548.5651 8.0039334 -2739.2844 190.71928 0.61 -2748.648 + 4000 1680.6126 58960.336 -2557.3154 8.0039334 -2744.7902 187.47487 0.6 -2753.7317 + 4100 1653.386 49607.708 -2575.5217 8.0039334 -2759.9594 184.43769 0.59 -2769.1633 + 4200 1628.5806 59615.376 -2562.3155 8.0039334 -2743.9861 181.67062 0.58 -2753.41 + 4300 1605.7638 57213.496 -2568.9241 8.0039334 -2748.0495 179.12536 0.57 -2758.0401 + 4400 1586.3539 57319.854 -2567.3656 8.0039334 -2744.3258 176.96016 0.56 -2754.0451 + 4500 1569.1762 51130.713 -2581.1996 8.0039334 -2756.2435 175.04397 0.55 -2767.0321 + 4600 1554.2279 51209.179 -2582.5415 8.0039334 -2755.918 173.37646 0.54 -2767.2662 + 4700 1541.2064 51391.73 -2584.9958 8.0039334 -2756.9197 171.92389 0.53 -2768.7565 + 4800 1529.3566 48164.466 -2590.6042 8.0039334 -2761.2062 170.60203 0.52 -2773.7243 + 4900 1517.8359 47270.397 -2595.8724 8.0039334 -2765.1893 169.31688 0.51 -2777.023 + 5000 1506.2157 46127.507 -2595.6789 8.0039334 -2763.6995 168.02063 0.5 -2777.2298 + 5100 1494.4675 45092.69 -2600.075 8.0039334 -2766.7851 166.7101 0.49 -2779.7181 + 5200 1480.2051 41485.775 -2607.0151 8.0039334 -2772.1342 165.1191 0.48 -2785.5483 + 5300 1465.1748 39849.775 -2612.3292 8.0039334 -2775.7716 163.44246 0.47 -2788.9493 + 5400 1447.9857 39330.898 -2615.0485 8.0039334 -2776.5735 161.52498 0.46 -2790.3996 + 5500 1427.9752 35674.647 -2624.4226 8.0039334 -2783.7154 159.29278 0.45 -2797.922 + 5600 1406.1138 33960.907 -2630.1163 8.0039334 -2786.9704 156.85411 0.44 -2801.8887 + 5700 1383.2834 39122.559 -2621.0951 8.0039334 -2775.4025 154.30734 0.43 -2790.3487 + 5800 1357.4953 34996.372 -2632.2759 8.0039334 -2783.7065 151.43064 0.42 -2798.8306 + 5900 1331.345 30380.598 -2640.9488 8.0039334 -2789.4624 148.51354 0.41 -2804.9873 + 6000 1304.5621 26910.214 -2652.7787 8.0039334 -2798.3046 145.52587 0.4 -2814.8623 + 6100 1278.2361 34273.343 -2642.9208 8.0039334 -2785.51 142.58917 0.39 -2801.9949 + 6200 1251.7315 29640.216 -2651.6543 8.0039334 -2791.2868 139.63253 0.38 -2807.555 + 6300 1226.3243 25302.881 -2660.8448 8.0039334 -2797.6431 136.79832 0.37 -2815.3412 + 6400 1202.0501 23978.218 -2668.9812 8.0039334 -2803.0717 134.09049 0.36 -2820.0498 + 6500 1178.656 25315.547 -2668.8266 8.0039334 -2800.3074 131.48085 0.35 -2817.3114 + 6600 1156.9362 22676.36 -2673.1256 8.0039334 -2802.1835 129.05797 0.34 -2819.6239 + 6700 1163.8953 22497.41 -2674.4164 8.0039334 -2804.2506 129.83427 0.33 -2821.3556 + 6800 1123.3867 18193.925 -2685.3997 8.0039334 -2810.7152 125.31548 0.32 -2830.2267 + 6900 1147.6164 19451.295 -2682.5112 8.0039334 -2810.5295 128.01834 0.31 -2830.1189 + 7000 1096.0437 16757.031 -2689.5567 8.0039334 -2811.822 122.26532 0.3 -2832.475 + 7100 1065.1232 11734.429 -2703.0123 8.0039334 -2821.8284 118.8161 0.29 -2841.8886 + 7200 1048.459 15187.579 -2697.6953 8.0039334 -2814.6525 116.95718 0.28 -2834.2885 + 7300 1067.9311 13697.721 -2698.4027 8.0039334 -2817.532 119.12932 0.27 -2838.9073 + 7400 998.4963 13330.213 -2707.915 8.0039334 -2819.2987 111.38377 0.26 -2840.2979 + 7500 1023.6105 7070.7227 -2712.3293 8.0039334 -2826.5146 114.18529 0.25 -2848.8706 + 7600 970.8256 7694.074 -2717.03 8.0039334 -2825.3271 108.29706 0.24 -2848.0154 + 7700 966.64807 2336.6244 -2727.1914 8.0039334 -2835.0224 107.83105 0.23 -2858.6468 + 7800 961.41051 564.58952 -2731.2826 8.0039334 -2838.5294 107.24679 0.22 -2862.4362 + 7900 880.63112 3569.7496 -2734.8614 8.0039334 -2833.0971 98.235727 0.21 -2857.2329 + 8000 904.4189 -883.03123 -2738.386 8.0039334 -2839.2753 100.88929 0.2 -2864.1778 + 8100 901.36129 -977.71992 -2739.797 8.0039334 -2840.3452 100.54821 0.19 -2864.7239 + 8200 834.53102 2535.4088 -2741.7378 8.0039334 -2834.831 93.093191 0.18 -2859.9007 + 8300 830.12299 -6039.909 -2755.2151 8.0039334 -2847.8166 92.601469 0.17 -2873.9613 + 8400 807.63492 -5386.015 -2755.7433 8.0039334 -2845.8362 90.092891 0.16 -2872.8494 + 8500 792.64609 -6372.4363 -2758.8515 8.0039334 -2847.2724 88.420864 0.15 -2874.8424 + 8600 758.58937 -11610.225 -2770.2002 8.0039334 -2854.822 84.621786 0.14 -2882.1413 + 8700 758.23023 -10178.597 -2769.9169 8.0039334 -2854.4986 84.581724 0.13 -2881.6683 + 8800 743.6008 -12026.268 -2776.0268 8.0039334 -2858.9766 82.949788 0.12 -2886.1335 + 8900 734.90694 -17472.954 -2783.3532 8.0039334 -2865.3332 81.979975 0.11 -2893.7244 + 9000 707.10989 -17006.609 -2785.8037 8.0039334 -2864.6829 78.879173 0.1 -2893.2911 + 9100 695.13334 -20282.325 -2793.2101 8.0039334 -2870.7532 77.543171 0.09 -2899.891 + 9200 642.26078 -20368.748 -2798.1031 8.0039334 -2869.7483 71.645157 0.08 -2899.7706 + 9300 641.79203 -21031.839 -2798.6175 8.0039334 -2870.2104 71.592866 0.07 -2901.0899 + 9400 640.84125 -20243.066 -2799.2216 8.0039334 -2870.7084 71.486807 0.06 -2901.4655 + 9500 609.96498 -23936.746 -2809.3896 8.0039334 -2877.4321 68.042512 0.05 -2908.2797 + 9600 574.1911 -26956.034 -2817.6178 8.0039334 -2881.6696 64.051882 0.04 -2912.6742 + 9700 567.47641 -28321.949 -2818.8602 8.0039334 -2882.163 63.302848 0.03 -2914.4223 + 9800 520.37758 -31000.6 -2828.0252 8.0039334 -2886.0741 58.048903 0.02 -2918.4921 + 9900 506.49368 -32223.744 -2830.1705 8.0039334 -2886.6706 56.500132 0.01 -2920.4741 + 10000 507.13597 -33930.476 -2834.3236 8.0039334 -2890.8954 56.571781 0 -2924.5423 +Loop time of 17.7279 on 2 procs for 10000 steps with 864 atoms + + +Total wall time: 0:00:17 diff --git a/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.1 b/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.1 new file mode 100644 index 0000000000..b71fad8cab --- /dev/null +++ b/examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.1 @@ -0,0 +1,172 @@ +LAMMPS (8 Feb 2023) +Processor partition = 1 + using 1 OpenMP thread(s) per MPI task +# Example for an alchemical transformation simulation of pure copper into an aluminium bronze + +# set up different names for two partitions +variable name world pure alloy + +# create pure copper system +units metal +lattice fcc 3.75 +Lattice spacing in x,y,z = 3.75 3.75 3.75 +region box block 0 6 0 6 0 6 +create_box 2 box +Created orthogonal box = (0 0 0) to (22.5 22.5 22.5) + 1 by 1 by 2 MPI processor grid + +timestep 0.002 +create_atoms 1 box +Created 864 atoms + using lattice units in orthogonal box = (0 0 0) to (22.5 22.5 22.5) + create_atoms CPU = 0.000 seconds +displace_atoms all random 0.3 0.3 0.3 57845645 +Displacing atoms ... +pair_style eam/alloy +pair_coeff * * AlCu.eam.alloy Cu Al +Reading eam/alloy potential file AlCu.eam.alloy with DATE: 2008-10-01 + +# replace 5% of copper with aluminium on the second partition only +if "${name} == alloy" then "set type 1 type/fraction 2 0.05 6745234" +set type 1 type/fraction 2 0.05 6745234 +Setting atom values ... + 37 settings made for type/fraction + +# define ramp variable to combine the two different partitions +if "${name} == pure" then "variable ramp equal ramp(1.0,0.0)" else "variable ramp equal ramp(0.0,1.0)" +variable ramp equal ramp(0.0,1.0) + +velocity all create 5000.0 6567345 +fix 1 all nvt temp 2500.0 500.0 0.002 +fix 2 all alchemy v_ramp + +compute pressure all pressure/alchemy 2 + +# only need to output a dump file from one partition +#if "${name} == alloy" then # "dump 1 all custom 100 ${name}.lammpstrj id type element x y z" # "dump_modify 1 element Cu Al" + +thermo_style custom step temp press etotal density pe ke f_2 f_2[3] +thermo_modify colname f_2 lambda colname f_2[3] EPot_mixed +thermo_modify press pressure + +thermo 100 +run 10000 post no +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 8.6825 + ghost atom cutoff = 8.6825 + binsize = 4.34125, bins = 6 6 6 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair eam/alloy, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 3.474 | 3.474 | 3.474 Mbytes + Step Temp Press TotEng Density PotEng KinEng lambda EPot_mixed + 0 5000 2135632.8 11134.401 7.8067107 10576.643 557.75753 0 10353.88 + 100 2502.4638 98461.445 -2410.2394 7.8067107 -2689.393 279.1536 0.01 -2673.3837 + 200 2471.2133 111478.55 -2390.1144 7.8067107 -2665.7819 275.66756 0.02 -2650.057 + 300 2442.869 109429.75 -2398.101 7.8067107 -2670.6067 272.50571 0.03 -2654.7075 + 400 2420.4514 97114.07 -2425.2882 7.8067107 -2695.2932 270.005 0.04 -2675.755 + 500 2403.3496 98285.746 -2421.1707 7.8067107 -2689.268 268.09727 0.05 -2673.6895 + 600 2388.6644 101125.13 -2418.9873 7.8067107 -2685.4464 266.45911 0.06 -2668.238 + 700 2372.4091 96916.491 -2428.7444 7.8067107 -2693.3902 264.64581 0.07 -2675.3394 + 800 2356.5786 101050.97 -2423.8648 7.8067107 -2686.7446 262.87989 0.08 -2669.6797 + 900 2337.3256 95290.542 -2439.2282 7.8067107 -2699.9604 260.73219 0.09 -2681.0381 + 1000 2314.6552 90348.199 -2452.6965 7.8067107 -2710.8998 258.20327 0.1 -2693.7969 + 1100 2288.9862 87089.44 -2457.7481 7.8067107 -2713.088 255.33986 0.11 -2696.9199 + 1200 2263.1553 89628.642 -2462.9772 7.8067107 -2715.4356 252.45838 0.12 -2698.4983 + 1300 2238.1696 92771.824 -2453.8691 7.8067107 -2703.5403 249.67118 0.13 -2687.6406 + 1400 2209.3333 83029.221 -2477.2462 7.8067107 -2723.7006 246.45446 0.14 -2706.638 + 1500 2181.9468 83386 -2476.3027 7.8067107 -2719.7022 243.39945 0.15 -2704.0915 + 1600 2156.3929 81801.467 -2483.9294 7.8067107 -2724.4783 240.54887 0.16 -2708.3304 + 1700 2133.221 84933.294 -2477.0832 7.8067107 -2715.0473 237.96401 0.17 -2700.142 + 1800 2110.4496 76700.779 -2495.3147 7.8067107 -2730.7386 235.42383 0.18 -2714.6978 + 1900 2091.8774 82568.985 -2488.7082 7.8067107 -2722.0603 233.35207 0.19 -2706.5643 + 2000 2074.8972 79024.983 -2496.7194 7.8067107 -2728.1773 231.45791 0.2 -2710.7 + 2100 2059.7023 72408.502 -2511.6849 7.8067107 -2741.4478 229.76289 0.21 -2725.1617 + 2200 2046.5649 81172.694 -2496.9796 7.8067107 -2725.277 228.2974 0.22 -2708.4618 + 2300 2033.9024 78029.129 -2500.8329 7.8067107 -2727.7178 226.88488 0.23 -2712.5894 + 2400 2022.0535 79691.571 -2499.8173 7.8067107 -2725.3804 225.56311 0.24 -2709.4434 + 2500 2008.7841 75015.359 -2512.1238 7.8067107 -2736.2067 224.08289 0.25 -2719.4159 + 2600 1995.175 74845.99 -2512.2807 7.8067107 -2734.8455 222.56477 0.26 -2718.1619 + 2700 1979.2542 76756.042 -2512.5095 7.8067107 -2733.2983 220.78879 0.27 -2716.0963 + 2800 1961.2518 70627.84 -2527.3853 7.8067107 -2746.1659 218.78059 0.28 -2730.7999 + 2900 1942.4706 70852.114 -2530.3101 7.8067107 -2746.9956 216.68552 0.29 -2731.5024 + 3000 1922.0973 75867.616 -2518.2668 7.8067107 -2732.6796 214.41284 0.3 -2718.4146 + 3100 1898.9538 68984.855 -2540.7552 7.8067107 -2752.5863 211.83116 0.31 -2738.0966 + 3200 1874.9107 73304.621 -2533.0578 7.8067107 -2742.2069 209.14911 0.32 -2727.1933 + 3300 1848.9824 69447.29 -2538.6006 7.8067107 -2744.8574 206.25677 0.33 -2730.768 + 3400 1821.4108 64973.868 -2552.8877 7.8067107 -2756.0688 203.18111 0.34 -2741.4173 + 3500 1795.4926 64688.566 -2555.5739 7.8067107 -2755.8638 200.2899 0.35 -2740.986 + 3600 1769.4009 56451.766 -2573.338 7.8067107 -2770.7174 197.37934 0.36 -2756.4626 + 3700 1745.5646 68268.147 -2554.1634 7.8067107 -2748.8838 194.72036 0.37 -2735.8618 + 3800 1721.4956 64187.24 -2566.0526 7.8067107 -2758.088 192.03542 0.38 -2742.67 + 3900 1699.517 61697.296 -2573.7098 7.8067107 -2763.2935 189.58368 0.39 -2748.648 + 4000 1680.2819 58960.336 -2579.7058 7.8067107 -2767.1438 187.43798 0.4 -2753.7317 + 4100 1662.2079 49607.708 -2596.9862 7.8067107 -2782.408 185.42179 0.41 -2769.1633 + 4200 1645.509 59615.376 -2582.865 7.8067107 -2766.424 183.55901 0.42 -2753.41 + 4300 1628.9336 57213.496 -2589.5734 7.8067107 -2771.2834 181.71 0.43 -2758.0401 + 4400 1614.2369 57319.854 -2586.3446 7.8067107 -2766.4151 180.07055 0.44 -2754.0451 + 4500 1598.6946 51130.713 -2601.8813 7.8067107 -2780.2181 178.33679 0.45 -2767.0321 + 4600 1583.2816 51209.179 -2603.9706 7.8067107 -2780.5881 176.61745 0.46 -2767.2662 + 4700 1567.3743 51391.73 -2607.2614 7.8067107 -2782.1044 174.84296 0.47 -2768.7565 + 4800 1550.9245 48164.466 -2614.2775 7.8067107 -2787.2854 173.00797 0.48 -2773.7243 + 4900 1532.697 47270.397 -2618.365 7.8067107 -2789.3397 170.97466 0.49 -2777.023 + 5000 1513.7604 46127.507 -2621.8979 7.8067107 -2790.7601 168.86225 0.5 -2777.2298 + 5100 1494.431 45092.69 -2625.438 7.8067107 -2792.144 166.70602 0.51 -2779.7181 + 5200 1472.5263 41485.775 -2633.6681 7.8067107 -2797.9306 164.26252 0.52 -2785.5483 + 5300 1450.7978 39849.775 -2638.7965 7.8067107 -2800.6352 161.83868 0.53 -2788.9493 + 5400 1428.063 39330.898 -2642.8749 7.8067107 -2802.1774 159.30258 0.54 -2790.3996 + 5500 1404.4139 35674.647 -2652.8811 7.8067107 -2809.5456 156.66449 0.55 -2797.922 + 5600 1380.8619 33960.907 -2659.5731 7.8067107 -2813.6103 154.03722 0.56 -2801.8887 + 5700 1357.9957 39122.559 -2650.1374 7.8067107 -2801.6239 151.48646 0.57 -2790.3487 + 5800 1334.8385 34996.372 -2660.8794 7.8067107 -2809.7826 148.90324 0.58 -2798.8306 + 5900 1312.6203 30380.598 -2669.3511 7.8067107 -2815.7758 146.42476 0.59 -2804.9873 + 6000 1290.7337 26910.214 -2681.9175 7.8067107 -2825.9008 143.98328 0.6 -2814.8623 + 6100 1269.6552 34273.343 -2670.9025 7.8067107 -2812.5345 141.63195 0.61 -2801.9949 + 6200 1249.4502 29640.216 -2678.1477 7.8067107 -2817.5258 139.37805 0.62 -2807.555 + 6300 1229.8377 25302.881 -2688.545 7.8067107 -2825.7352 137.19025 0.63 -2815.3412 + 6400 1211.0142 23978.218 -2694.5096 7.8067107 -2829.6 135.09045 0.64 -2820.0498 + 6500 1192.3094 25315.547 -2693.4635 7.8067107 -2826.4674 133.0039 0.65 -2817.3114 + 6600 1174.741 22676.36 -2697.5642 7.8067107 -2828.6083 131.04413 0.66 -2819.6239 + 6700 1117.0829 22497.41 -2705.1682 7.8067107 -2829.7804 124.61228 0.67 -2821.3556 + 6800 1125.4062 18193.925 -2713.8679 7.8067107 -2839.4086 125.54075 0.68 -2830.2267 + 6900 1142.0728 19451.295 -2711.5199 7.8067107 -2838.9199 127.39994 0.69 -2830.1189 + 7000 1148.9469 16757.031 -2713.1595 7.8067107 -2841.3263 128.16676 0.7 -2832.475 + 7100 1044.2332 11734.429 -2733.5965 7.8067107 -2850.0822 116.48578 0.71 -2841.8886 + 7200 1052.3967 15187.579 -2724.5284 7.8067107 -2841.9248 117.39644 0.72 -2834.2885 + 7300 1045.9123 13697.721 -2730.1401 7.8067107 -2846.8132 116.67309 0.73 -2838.9073 + 7400 1082.9806 13330.213 -2726.8679 7.8067107 -2847.676 120.80812 0.74 -2840.2979 + 7500 980.30736 7070.7227 -2746.9678 7.8067107 -2856.3226 109.35476 0.75 -2848.8706 + 7600 982.86602 7694.074 -2745.5399 7.8067107 -2855.1801 109.64018 0.76 -2848.0154 + 7700 958.27863 2336.6244 -2758.806 7.8067107 -2865.7035 106.89742 0.77 -2858.6468 + 7800 991.25495 564.58952 -2758.6032 7.8067107 -2869.1791 110.57598 0.78 -2862.4362 + 7900 901.65217 3569.7496 -2763.0681 7.8067107 -2863.6488 100.58066 0.79 -2857.2329 + 8000 871.00782 -883.03123 -2773.2412 7.8067107 -2870.4035 97.162233 0.8 -2864.1778 + 8100 873.73598 -977.71992 -2772.9758 7.8067107 -2870.4424 97.466563 0.81 -2864.7239 + 8200 859.27838 2535.4088 -2769.5501 7.8067107 -2865.4039 95.853797 0.82 -2859.9007 + 8300 874.88664 -6039.909 -2781.7213 7.8067107 -2879.3163 97.594921 0.83 -2873.9613 + 8400 810.33098 -5386.015 -2787.6011 7.8067107 -2877.9947 90.393641 0.84 -2872.8494 + 8500 812.50682 -6372.4363 -2789.0713 7.8067107 -2879.7077 90.636359 0.85 -2874.8424 + 8600 792.27145 -11610.225 -2798.2095 7.8067107 -2886.5886 88.379073 0.86 -2882.1413 + 8700 737.92522 -10178.597 -2803.4114 7.8067107 -2885.7281 82.316669 0.87 -2881.6683 + 8800 744.31065 -12026.268 -2806.8078 7.8067107 -2889.8367 83.028973 0.88 -2886.1335 + 8900 709.31385 -17472.954 -2818.1084 7.8067107 -2897.2334 79.125027 0.89 -2893.7244 + 9000 704.20763 -17006.609 -2817.9144 7.8067107 -2896.4698 78.555422 0.9 -2893.2911 + 9100 697.72956 -20282.325 -2824.94 7.8067107 -2902.7728 77.832782 0.91 -2899.891 + 9200 696.15397 -20368.748 -2824.7243 7.8067107 -2902.3813 77.657023 0.92 -2899.7706 + 9300 627.98752 -21031.839 -2833.3613 7.8067107 -2903.4142 70.052953 0.93 -2901.0899 + 9400 603.4933 -20243.066 -2836.1082 7.8067107 -2903.4288 67.320586 0.94 -2901.4655 + 9500 599.69375 -23936.746 -2843.0066 7.8067107 -2909.9033 66.89674 0.95 -2908.2797 + 9600 587.42611 -26956.034 -2848.4377 7.8067107 -2913.966 65.528267 0.96 -2912.6742 + 9700 567.47718 -28321.949 -2852.1171 7.8067107 -2915.42 63.302933 0.97 -2914.4223 + 9800 524.81105 -31000.6 -2860.6103 7.8067107 -2919.1537 58.543462 0.98 -2918.4921 + 9900 507.60008 -32223.744 -2864.192 7.8067107 -2920.8155 56.623552 0.99 -2920.4741 + 10000 521.28866 -33930.476 -2866.3918 7.8067107 -2924.5423 58.150534 1 -2924.5423 +Loop time of 17.7277 on 2 procs for 10000 steps with 864 atoms + + +Total wall time: 0:00:17 diff --git a/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4 b/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4 new file mode 100644 index 0000000000..18375012d0 --- /dev/null +++ b/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4 @@ -0,0 +1,13 @@ +LAMMPS (8 Feb 2023) +Running on 2 partitions of processors +Starting alchemical transformation at 0% + Alchemical transformation progress: 10% + Alchemical transformation progress: 20% + Alchemical transformation progress: 30% + Alchemical transformation progress: 40% + Alchemical transformation progress: 50% + Alchemical transformation progress: 60% + Alchemical transformation progress: 70% + Alchemical transformation progress: 80% + Alchemical transformation progress: 90% + Alchemical transformation progress: 100% diff --git a/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.0 b/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.0 new file mode 100644 index 0000000000..49e6c6536d --- /dev/null +++ b/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.0 @@ -0,0 +1,412 @@ +LAMMPS (8 Feb 2023) +Processor partition = 0 + using 1 OpenMP thread(s) per MPI task +# Example for an alchemical transformation of two water molecules into a hydronium and hydroxyl ion +# WARNING: This input is intended for demonstrating the method only, +# the force field parameters are mostly made up and NOT suitable for production simulations. + +# set up different names for two partitions +variable name world twowater twoions + +units real +atom_style full +atom_modify map array +region box block -5 5 -5 5 -5 5 +boundary p p p +create_box 2 box bond/types 2 angle/types 2 extra/bond/per/atom 3 extra/angle/per/atom 3 extra/special/per/atom 3 +Created orthogonal box = (-5 -5 -5) to (5 5 5) + 1 by 1 by 2 MPI processor grid + +mass 1 15.9994 +mass 2 1.008 + +pair_style lj/cut/coul/cut 10.0 +pair_coeff 1 1 0.1553 3.166 +pair_coeff 1 2 0.0 1.0 +pair_coeff 2 2 0.0 1.0 + +bond_style harmonic +bond_coeff * 1000.0 1.0 + +angle_style harmonic +angle_coeff * 100.0 109.47 + +molecule water h2o.mol +Read molecule template water: + 1 molecules + 0 fragments + 3 atoms with max type 2 + 2 bonds with max type 1 + 1 angles with max type 1 + 0 dihedrals with max type 0 + 0 impropers with max type 0 + +# create the two molecules we want to transform ... +create_atoms 0 single -2.0 0.0 0.0 mol water 453624 +Created 3 atoms + using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5) + create_atoms CPU = 0.003 seconds +create_atoms 0 single 2.0 0.0 0.0 mol water 767353 +Created 3 atoms + using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5) + create_atoms CPU = 0.000 seconds + +# ... and put them in a group +group transform id 1:6 +6 atoms in group transform + +# now fill the rest of the box with more water +create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33 +Created 96 atoms + using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5) + create_atoms CPU = 0.000 seconds + +# change topology and settings for the two states +# we cannot simply create a different topology directly or +# load a different data file because the order and position +# of all atoms must be maintained across both replica + +# we first have to remove all topology data in the transform group +delete_bonds transform bond 1 +System init for delete_bonds ... +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 12 + ghost atom cutoff = 12 + binsize = 6, bins = 2 2 2 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair lj/cut/coul/cut, perpetual + attributes: half, newton on + pair build: half/bin/newton + stencil: half/bin/3d + bin: standard +Deleting bonds ... + 68 total bonds, 64 turned on, 4 turned off + 34 total angles, 34 turned on, 0 turned off + 0 total dihedrals, 0 turned on, 0 turned off + 0 total impropers, 0 turned on, 0 turned off +delete_bonds transform angle 1 remove +System init for delete_bonds ... +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +Deleting bonds ... + 64 total bonds, 64 turned on, 0 turned off + 32 total angles, 32 turned on, 0 turned off + 0 total dihedrals, 0 turned on, 0 turned off + 0 total impropers, 0 turned on, 0 turned off + +# then generate different topologies for the two partitions. select by name. +if "${name} == twowater" then "create_bonds single/bond 2 1 2" "create_bonds single/bond 2 1 3" "create_bonds single/bond 2 4 5" "create_bonds single/bond 2 4 6" "create_bonds single/angle 2 2 1 3" "create_bonds single/angle 2 5 4 6" else "create_bonds single/bond 2 1 2" "create_bonds single/bond 2 3 4" "create_bonds single/bond 2 4 5" "create_bonds single/bond 2 4 6" "create_bonds single/angle 2 3 4 5" "create_bonds single/angle 2 5 4 6" "create_bonds single/angle 2 3 4 6" "set atom 1 charge -1.1354" "set atom 2 charge 0.1354" "set atom 3 charge 0.56775" "set atom 4 charge -0.70305" "set atom 5*6 charge 0.56775" +create_bonds single/bond 2 1 2 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/bond 2 1 3 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/bond 2 4 5 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/bond 2 4 6 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/angle 2 2 1 3 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/angle 2 5 4 6 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds + +velocity all create 300.0 5463576 +timestep 0.2 + +# define ramp variable to combine the two different partitions +if "${name} == twowater" then "variable ramp equal ramp(1.0,0.0)" else "variable ramp equal ramp(0.0,1.0)" +variable ramp equal ramp(1.0,0.0) + +# since the trajectory and forces are kept identical through fix alchemy, +# we can do fix npt simulations, but we must use the "mixed" pressure + +fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0 +fix transform all alchemy v_ramp +compute pressure all pressure/alchemy transform +fix_modify integrate press pressure + +# only need to output a dump file from one partition +# if "${name} == twowater" then # "dump 1 all atom 100 ${name}.lammpstrj" # "dump_modify 1 sort id" + +thermo_style custom step temp press etotal density pe ke f_transform f_transform[3] +thermo_modify colname f_transform lambda colname f_transform[3] EPot_mixed +thermo_modify press pressure + +thermo 100 +run 20000 +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +Per MPI rank memory allocation (min/avg/max) = 7.535 | 7.535 | 7.535 Mbytes + Step Temp Press TotEng Density PotEng KinEng lambda EPot_mixed + 0 300 184631.17 -802.19609 1.0171193 -892.51464 90.318545 1 -892.51464 + 100 293.4335 -513.97927 189.09028 0.3029684 100.74865 88.341623 0.995 207.29243 + 200 320.81838 1513.3 111.09325 0.64509407 14.507087 96.586165 0.99 46.402434 + 300 261.53825 -3808.6942 -538.67351 0.88817871 -617.41269 78.739179 0.985 -614.76481 + 400 308.32444 1375.0327 -804.08042 0.97694485 -896.90513 92.824715 0.98 -897.24386 + 500 292.60623 -12843.384 -1141.5985 1.0140814 -1229.691 88.092563 0.975 -1221.948 + 600 326.16112 91608.645 2882.4855 1.0263567 2784.2908 98.194659 0.97 2786.0743 + 700 317.22811 -60420.408 -1517.0601 0.9919143 -1612.5654 95.505272 0.965 -1616.7677 + 800 309.94673 -56793.931 -1549.4612 1.0024958 -1642.7743 93.313126 0.96 -1640.8606 + 900 272.80718 -5076.2282 -1266.6775 1.0141078 -1348.8093 82.131824 0.955 -1357.1812 + 1000 346.52557 -16919.192 -1406.675 1.0170322 -1511.0006 104.32562 0.95 -1514.1398 + 1100 300.62984 -51826.972 -2469.2821 1.0063592 -2559.7903 90.508164 0.945 -2565.9668 + 1200 293.72767 45222.698 2039.1938 1.0272336 1950.7636 88.430186 0.94 1951.1683 + 1300 289.4971 94600.203 2184.5096 1.0298091 2097.3531 87.156524 0.935 2088.1011 + 1400 328.98666 35468.53 1816.723 1.0198467 1717.6777 99.045321 0.93 1716.9268 + 1500 392.9407 41220.21 1820.7312 1.0194044 1702.4317 118.29944 0.925 1690.6416 + 1600 284.90093 -37995.141 -2439.9256 1.0068497 -2525.6984 85.772791 0.92 -2529.5874 + 1700 267.67454 25275.255 1612.5912 1.0313556 1532.0046 80.586582 0.915 1520.9695 + 1800 247.16495 -68222.994 -2699.1903 1.0018048 -2773.6022 74.41193 0.91 -2788.9545 + 1900 328.41849 -51979.665 -2339.7432 1.0154952 -2438.6175 98.874268 0.905 -2461.3422 + 2000 313.8715 21726.003 1473.9594 1.018831 1379.4647 94.494725 0.9 1363.0739 + 2100 296.07341 66130.11 1379.1597 1.0224882 1290.0233 89.136398 0.895 1272.1884 + 2200 282.67889 51746.464 2376.377 1.0362683 2291.2732 85.103821 0.89 2276.1488 + 2300 281.04207 -71937.759 -2791.0867 1.0111968 -2875.6977 84.611035 0.885 -2892.2526 + 2400 321.00806 -51658.758 -2507.6376 1.0168161 -2604.2809 96.643269 0.88 -2622.9732 + 2500 288.95906 49645.753 2361.024 1.0355648 2274.0295 86.99454 0.875 2275.0348 + 2600 303.0064 -48470.216 -1973.9561 1.000439 -2065.1798 91.223657 0.87 -2082.5895 + 2700 361.7083 54973.584 1929.0772 1.0223169 1820.1806 108.89656 0.865 1817.1972 + 2800 289.27827 33564.017 1783.5154 1.0257088 1696.4248 87.090642 0.86 1682.9162 + 2900 283.94365 -47190.274 -2286.9529 1.0059224 -2372.4375 85.484591 0.855 -2395.9575 + 3000 337.10996 61988.686 2289.363 1.0258 2187.8721 101.49094 0.85 2168.1085 + 3100 305.28249 -50123.109 -1913.7375 1.0143567 -2005.6464 91.9089 0.845 -2040.3453 + 3200 270.19471 -18114.265 -1961.3048 1.0132413 -2042.6502 81.345309 0.84 -2076.9595 + 3300 263.39996 -109301.89 -2351.6647 0.99956867 -2430.9643 79.299671 0.835 -2471.0269 + 3400 249.01722 -37585.949 -2423.3849 1.0080932 -2498.3544 74.969576 0.83 -2537.734 + 3500 267.88868 65528.865 2164.3882 1.0390872 2083.7372 80.651054 0.825 2061.6421 + 3600 351.46458 -67583.574 -2861.4561 1.0059621 -2967.2687 105.81256 0.82 -3004.5028 + 3700 332.82565 -61204.843 -2275.2765 1.0053485 -2375.4776 100.20109 0.815 -2415.4435 + 3800 298.0279 61123.085 1951.1943 1.020672 1861.4695 89.724819 0.81 1831.9734 + 3900 317.04715 43210.045 1848.3875 1.0287838 1752.9367 95.450792 0.805 1721.1961 + 4000 261.75174 -111920.72 -3274.707 0.99457106 -3353.5105 78.803455 0.8 -3394.3528 + 4100 338.34245 89763.083 2115.4809 1.0318272 2013.6189 101.86199 0.795 1973.419 + 4200 288.52623 19989.792 1758.0617 1.0204119 1671.1975 86.86423 0.79 1633.1439 + 4300 322.05593 -33791.403 -2107.1643 1.0155065 -2204.1231 96.958741 0.785 -2246.4353 + 4400 280.59054 -27839.599 -1534.3265 1.0123506 -1618.8016 84.475096 0.78 -1672.1425 + 4500 296.08336 75206.112 2789.7494 1.0329133 2700.61 89.139393 0.775 2670.4908 + 4600 288.02094 -67445.101 -2983.5242 0.99579107 -3070.2363 86.712106 0.77 -3132.6352 + 4700 332.21641 -83964.647 -2746.0748 1.0053915 -2846.0925 100.01768 0.765 -2916.6652 + 4800 305.52557 -93111.682 -2354.6562 0.99824941 -2446.6383 91.982082 0.76 -2514.039 + 4900 272.82262 -53664.198 -2790.1342 0.99527907 -2872.2707 82.136472 0.755 -2947.2062 + 5000 288.96388 39846.089 1975.7895 1.0371205 1888.7935 86.99599 0.75 1833.3716 + 5100 293.25543 -71819.425 -2916.8817 0.99452016 -3005.1697 88.288011 0.745 -3069.7618 + 5200 296.16486 -64484.821 -2416.0381 1.0153125 -2505.2021 89.163931 0.74 -2589.6878 + 5300 342.07246 60255.985 1749.8984 1.0306374 1646.9134 102.98496 0.735 1608.5673 + 5400 295.53315 61421.189 1620.6432 1.0334762 1531.6695 88.973746 0.73 1492.0969 + 5500 285.1755 24845.113 1144.9009 1.0335035 1059.0454 85.855455 0.725 1014.0263 + 5600 276.79246 66244.072 1523.3405 1.0251134 1440.0089 83.331641 0.72 1390.9078 + 5700 288.43485 14629.129 1208.7278 1.02997 1121.8911 86.836719 0.715 1074.3779 + 5800 276.64233 59262.97 914.26981 1.0332305 830.98336 83.286443 0.71 763.27415 + 5900 358.25384 -111203.23 -3370.7104 1.006331 -3478.567 107.85655 0.705 -3561.6283 + 6000 327.82233 -62936.335 -1881.2426 1.0106476 -1979.9373 98.694786 0.7 -2062.107 + 6100 309.86553 -59110.381 -2698.4147 1.0076375 -2791.7034 93.28868 0.695 -2870.9566 + 6200 294.95752 -58011.845 -1819.201 1.0161986 -1908.0014 88.800448 0.69 -1971.2586 + 6300 307.27252 59016.127 1186.0709 1.0225343 1093.5629 92.508022 0.685 1032.9616 + 6400 292.21611 -118857.83 -2992.9227 1.0031304 -3080.8978 87.975114 0.68 -3175.9068 + 6500 332.21108 97132.547 2290.2726 1.0349752 2190.2565 100.01607 0.675 2121.4674 + 6600 321.93985 -50017.836 -2484.5919 1.0076385 -2581.5157 96.923796 0.67 -2673.6819 + 6700 299.07231 -98136.241 -3013.1584 1.0106261 -3103.1976 90.039252 0.665 -3194.0694 + 6800 285.33731 28487.407 455.68104 1.017186 369.77687 85.904169 0.66 294.02418 + 6900 296.41387 50054.533 2028.1018 1.0445934 1938.8629 89.238897 0.655 1851.1397 + 7000 318.34109 -71579.897 -3114.0324 1.007304 -3209.8728 95.840345 0.65 -3305.5715 + 7100 317.2192 12370.315 1139.7687 1.0177238 1044.2661 95.502589 0.645 961.68146 + 7200 289.91091 43944.333 1213.4685 1.0302599 1126.1874 87.281105 0.64 1039.3329 + 7300 324.81227 42859.173 1288.5297 1.0300698 1190.7411 97.788572 0.635 1083.5764 + 7400 305.72634 76516.72 2404.9462 1.0580673 2312.9036 92.042527 0.63 2218.3674 + 7500 283.77479 -87225.247 -3401.4354 0.9936609 -3486.8692 85.433754 0.625 -3581.9003 + 7600 318.72048 16084.995 1848.3396 1.0296163 1752.385 95.954566 0.62 1653.5125 + 7700 307.36895 74007.024 1754.8928 1.0363553 1662.3558 92.537054 0.615 1572.8568 + 7800 297.44704 28693.429 1261.5958 1.0206402 1172.0459 89.549946 0.61 1093.5903 + 7900 288.47173 -94644.76 -2455.4369 1.0095244 -2542.2847 86.847823 0.605 -2648.2974 + 8000 280.5001 63255.987 1584.7287 1.020587 1500.2808 84.44787 0.6 1418.7144 + 8100 320.23351 68383.584 2300.142 1.0419027 2203.7319 96.410082 0.595 2111.2616 + 8200 354.86765 50012.316 1350.3752 1.0292074 1243.5381 106.8371 0.59 1149.0631 + 8300 278.8859 24200.145 772.82478 1.0232468 688.86288 83.961894 0.585 569.83331 + 8400 295.95754 -12010.401 1078.6687 1.0200526 989.56716 89.101514 0.58 872.63254 + 8500 289.57525 -64114.236 -3017.5291 1.0062613 -3104.7091 87.18005 0.575 -3227.8592 + 8600 329.23192 15948.893 1427.4791 1.0271511 1328.3599 99.119161 0.57 1211.3868 + 8700 327.45998 39407.811 1697.908 1.0343756 1599.3223 98.585695 0.565 1480.9482 + 8800 296.99753 18017.036 871.11717 1.0178295 781.70256 89.414616 0.56 657.09357 + 8900 289.10016 79527.35 2536.8041 1.0462904 2449.7671 87.037018 0.555 2338.4742 + 9000 260.53308 -51430.717 -2680.1664 1.0040128 -2758.603 78.436561 0.55 -2894.7913 + 9100 314.17867 40043.249 2043.6144 1.0361707 1949.0272 94.5872 0.545 1836.5399 + 9200 362.21236 -55153.163 -2232.4614 1.0102652 -2341.5097 109.04831 0.54 -2497.8765 + 9300 266.62049 -72193.959 -2980.9254 1.0063615 -3061.1946 80.269248 0.535 -3225.2604 + 9400 259.59316 7674.9295 1378.2484 1.0326879 1300.0948 78.153588 0.53 1171.7119 + 9500 278.72962 -33127.835 -2581.8908 1.0163553 -2665.8057 83.914845 0.525 -2805.1711 + 9600 299.09039 24570.664 1839.3275 1.025317 1749.2828 90.044695 0.52 1631.8392 + 9700 336.65655 37767.059 1415.9074 1.0245842 1314.5529 101.35443 0.515 1218.4022 + 9800 279.68112 -55007.939 -2548.4947 1.0129837 -2632.696 84.201305 0.51 -2775.6909 + 9900 300.55694 14964.424 1332.287 1.0176315 1241.8008 90.486218 0.505 1127.7276 + 10000 276.77193 -75638.817 -3100.7483 1.0019055 -3184.0738 83.325459 0.5 -3341.663 + 10100 316.26926 61548.301 2406.8933 1.0427988 2311.6767 95.216598 0.495 2206.7159 + 10200 297.28171 56814.839 2073.2348 1.0270036 1983.7347 89.500171 0.49 1859.0176 + 10300 276.84774 83572.576 1920.2202 1.0322243 1836.8719 83.348283 0.485 1704.2233 + 10400 263.28687 -16549.994 1144.9186 1.0177296 1065.6529 79.265623 0.48 940.36339 + 10500 332.79252 47379.311 1538.8168 1.0223403 1438.6257 100.19112 0.475 1301.3762 + 10600 329.01105 39299.248 1330.9914 1.0230094 1231.9387 99.052663 0.47 1077.5639 + 10700 339.21335 -4217.0313 981.9193 1.0235915 879.79511 102.12419 0.465 741.87272 + 10800 297.11192 -80013.424 -3682.9794 1.0009566 -3772.4285 89.449054 0.46 -3958.8283 + 10900 299.63096 47783.864 1833.1832 1.0435283 1742.9757 90.207441 0.455 1620.9328 + 11000 363.58243 67370.65 1769.9236 1.0450565 1660.4628 109.46079 0.45 1543.803 + 11100 306.04664 50246.004 1323.7676 1.0341876 1231.6287 92.138957 0.445 1063.5345 + 11200 255.68488 43750.967 2205.2392 1.0530607 2128.2623 76.976955 0.44 1986.9636 + 11300 271.71957 -40289.299 -3083.4093 1.0168497 -3165.2137 81.804388 0.435 -3344.1558 + 11400 291.74007 -115249.24 -3319.2153 1.0067921 -3407.0471 87.831795 0.43 -3569.0256 + 11500 270.22077 28235.119 596.78348 1.0195886 515.43032 81.353155 0.425 373.78255 + 11600 297.08506 29790.701 1031.3514 1.0317561 941.91046 89.440969 0.42 803.92835 + 11700 295.33119 -101140.58 -3068.1166 1.0107252 -3157.0295 88.912943 0.415 -3348.5652 + 11800 308.77287 -120312.75 -3553.1266 0.99681611 -3646.0863 92.95972 0.41 -3835.4663 + 11900 293.86337 82287.969 1847.0953 1.0331446 1758.6242 88.47104 0.405 1625.4846 + 12000 279.9866 27919.689 2249.3175 1.0384909 2165.0242 84.293273 0.4 2028.6311 + 12100 292.19439 47035.018 1220.3593 1.0210948 1132.3908 87.968573 0.395 976.68113 + 12200 317.78325 46886.483 2039.3721 1.0310567 1943.6997 95.672401 0.39 1763.6359 + 12300 357.44038 51492.35 1898.3396 1.0414835 1790.7279 107.61165 0.385 1638.1905 + 12400 288.16042 -50528.478 -2837.5599 1.0130798 -2924.314 86.7541 0.38 -3115.1383 + 12500 300.37356 -114034.29 -3091.586 1.0130272 -3182.017 90.431009 0.375 -3364.511 + 12600 265.4402 73960.28 1297.3012 1.0343798 1217.3873 79.913907 0.37 1094.544 + 12700 287.07994 22722.337 1711.9055 1.0268893 1625.4767 86.428809 0.365 1490.443 + 12800 315.39745 46072.191 926.21362 1.0205906 831.25949 94.954129 0.36 662.83753 + 12900 307.48549 -120472.73 -3181.5418 1.005689 -3274.1139 92.57214 0.355 -3469.7834 + 13000 291.00577 -59154.878 -2652.0402 1.012821 -2739.6509 87.610727 0.35 -2966.2448 + 13100 318.37718 86703.775 2087.145 1.036209 1991.2938 95.851211 0.345 1819.7078 + 13200 280.41509 -107033.79 -2683.9682 1.0135499 -2768.3905 84.422277 0.34 -2953.8355 + 13300 275.01919 112086.85 2632.1169 1.0437945 2549.3192 82.797775 0.335 2394.6707 + 13400 299.46547 -111959.62 -3475.4354 0.98732074 -3565.593 90.157617 0.33 -3792.6798 + 13500 275.65631 17625.49 1762.6531 1.0255473 1679.6636 82.989588 0.325 1521.2745 + 13600 302.63525 122466.63 2853.8095 1.0438137 2762.6976 91.111917 0.32 2611.7188 + 13700 302.21323 -94399.663 -3003.7573 1.0068914 -3094.7422 90.984863 0.315 -3328.4539 + 13800 305.70483 30039.825 1792.3286 1.0225478 1700.2926 92.036052 0.31 1502.8476 + 13900 273.15564 -27822.597 -2346.6702 1.0131002 -2428.9069 82.236732 0.305 -2683.7804 + 14000 257.14635 -1821.0012 1378.2522 1.0172653 1300.8352 77.416947 0.3 1067.8821 + 14100 323.74174 -76139.363 -2747.4152 1.0139806 -2844.8814 97.466276 0.295 -3052.598 + 14200 274.99913 54440.943 983.27131 1.0186356 900.47957 82.791736 0.29 767.64274 + 14300 291.8774 1281.7286 1258.7881 1.0214356 1170.915 87.873139 0.285 1042.8585 + 14400 308.97245 -47961.551 -2680.7892 1.0155221 -2773.809 93.019806 0.28 -3007.3204 + 14500 311.05309 86026.739 2115.2413 1.0311818 2021.5951 93.646208 0.275 1830.5463 + 14600 269.27432 -12972.747 971.48792 1.0172239 890.41971 81.068214 0.27 711.81484 + 14700 262.20739 -77722.459 -3215.6796 1.0029421 -3294.6202 78.940634 0.265 -3537.1929 + 14800 300.33592 -107332.8 -2798.6664 1.0126785 -2889.086 90.419677 0.26 -3125.8439 + 14900 307.0173 -73204.94 -2652.7275 1.016875 -2745.1587 92.431185 0.255 -3016.0143 + 15000 334.08711 64440.709 1615.3385 1.0369481 1514.7576 100.58087 0.25 1323.9179 + 15100 293.46703 35179.395 2020.1996 1.0331703 1931.8479 88.351717 0.245 1768.9681 + 15200 307.80546 54006.291 1111.7619 1.0246134 1019.0935 92.668471 0.24 822.22754 + 15300 336.59792 33832.907 1282.8271 1.0242407 1181.4903 101.33678 0.235 956.25285 + 15400 285.46526 -99629.94 -2336.7535 1.0041261 -2422.6962 85.942688 0.23 -2730.1887 + 15500 269.15938 122546.45 3012.0792 1.0318941 2931.0456 81.033612 0.225 2695.4331 + 15600 294.50818 -104733.19 -1914.0242 1.0017241 -2002.6894 88.665166 0.22 -2281.5773 + 15700 321.05125 -40438.501 -1721.9536 1.0125007 -1818.6099 96.656271 0.215 -2084.0893 + 15800 303.22185 -42224.408 -2145.481 1.0090867 -2236.7695 91.288521 0.21 -2506.3355 + 15900 281.65025 78652.063 2791.9418 1.0460521 2707.1477 84.794136 0.205 2500.7128 + 16000 286.5056 63174.412 1725.0409 1.0185277 1638.785 86.255897 0.2 1442.2427 + 16100 271.74574 -75051.584 -2101.5778 1.0091727 -2183.3901 81.812265 0.195 -2447.1967 + 16200 266.40007 -64120.09 -2223.3693 1.0064388 -2303.5722 80.202888 0.19 -2525.4243 + 16300 330.67137 46204.124 1837.3168 1.0271037 1737.7643 99.552523 0.185 1529.8186 + 16400 310.94573 18510.676 1531.823 1.0188695 1438.2091 93.613887 0.18 1195.3041 + 16500 270.29191 79005.409 1989.9392 1.0354088 1908.5646 81.374574 0.175 1678.7342 + 16600 334.24108 17143.221 1644.1483 1.0230561 1543.521 100.62723 0.17 1379.3215 + 16700 292.28721 -78618.031 -2768.3997 0.99940789 -2856.3962 87.996518 0.165 -3087.379 + 16800 346.09157 -69398.937 -2647.1375 1.0153026 -2751.3325 104.19496 0.16 -3011.3476 + 16900 284.37042 39314.526 2287.7061 1.0343186 2202.093 85.613074 0.155 2045.8316 + 17000 281.69389 86032.163 1817.9125 1.0237583 1733.1052 84.807275 0.15 1528.8391 + 17100 297.2629 55128.393 2617.4156 1.0419182 2527.9211 89.494508 0.145 2369.0347 + 17200 328.81472 47063.86 1880.9768 1.0364627 1781.9832 98.993555 0.14 1526.7115 + 17300 313.41363 -46685.063 -2310.803 1.008117 -2405.1599 94.356878 0.135 -2685.8034 + 17400 288.17114 68721.523 2771.6788 1.0429105 2684.9215 86.757326 0.13 2490.5845 + 17500 260.01702 -23324.539 -1436.6738 1.0100415 -1514.955 78.281195 0.125 -1812.4813 + 17600 297.6633 72795.463 2635.741 1.0362827 2546.126 89.615054 0.12 2336.0152 + 17700 301.02937 93833.118 2690.518 1.0276466 2599.8896 90.628448 0.115 2364.8751 + 17800 238.92837 -78317.212 -2175.2589 1.0084077 -2247.1911 71.932208 0.11 -2579.4429 + 17900 299.27724 -65278.891 -2589.1173 1.0034773 -2679.2182 90.10095 0.105 -2987.498 + 18000 315.15753 -72437.481 -2338.6068 1.0136139 -2433.4887 94.881898 0.1 -2718.125 + 18100 268.39411 -42194.209 -2285.3824 1.0108104 -2366.1856 80.803219 0.095 -2656.3853 + 18200 284.91941 40905.162 2045.1299 1.0266185 1959.3515 85.778354 0.09 1723.6803 + 18300 329.41629 -56526.814 -2104.8386 1.0137784 -2204.0133 99.174667 0.085 -2540.6829 + 18400 293.16732 -22085.551 -1652.5444 1.0163557 -1740.8059 88.261484 0.08 -2052.7864 + 18500 276.28141 29367.359 2314.5497 1.0270087 2231.3719 83.177784 0.075 1971.9405 + 18600 331.18031 94669.124 2466.3928 1.0225403 2366.6871 99.705745 0.07 2081.8675 + 18700 265.4491 49872.536 2571.9864 1.0255308 2492.0698 79.916587 0.065 2260.4261 + 18800 314.51042 60758.629 2560.3258 1.0300818 2465.6387 94.687078 0.06 2219.6547 + 18900 317.70309 66345.835 1768.2041 1.0248229 1672.5558 95.648268 0.055 1410.3087 + 19000 291.39831 -96336.931 -2448.6251 1.0100498 -2536.354 87.728905 0.05 -2805.9396 + 19100 303.37353 -21421.495 -2012.0929 0.99771284 -2103.4271 91.334186 0.045 -2436.2173 + 19200 293.92948 36383.881 2042.0163 1.0251299 1953.5254 88.490944 0.04 1633.579 + 19300 287.24409 -46691.027 -2399.698 1.0112067 -2486.1762 86.478228 0.035 -2811.3506 + 19400 333.42205 21152.63 1524.8699 1.0189912 1424.4893 100.38065 0.03 1118.2463 + 19500 309.17989 60445.761 2298.6966 1.0338543 2205.6144 93.082259 0.025 1928.6181 + 19600 255.92246 -32795.547 -1724.7933 1.0157044 -1801.8418 77.04848 0.02 -2213.9484 + 19700 342.62332 33538.746 2439.138 1.0225905 2335.9872 103.1508 0.015 2081.9185 + 19800 295.68946 81474.342 2308.3249 1.0240743 2219.3041 89.020804 0.01 1910.5545 + 19900 306.4947 43488.052 2330.4958 1.0222707 2238.222 92.27385 0.005 1957.4368 + 20000 313.31679 -25133.284 -1161.6979 1.0163289 -1256.0256 94.327722 0 -1649.7551 +Loop time of 8.92653 on 2 procs for 20000 steps with 102 atoms + +Performance: 38.716 ns/day, 0.620 hours/ns, 2240.513 timesteps/s, 228.532 katom-step/s +96.3% CPU use with 2 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.7448 | 5.9908 | 7.2367 | 50.9 | 67.11 +Bond | 0.036126 | 0.03641 | 0.036695 | 0.1 | 0.41 +Neigh | 0.082247 | 0.082294 | 0.082342 | 0.0 | 0.92 +Comm | 0.68467 | 1.9303 | 3.1759 | 89.7 | 21.62 +Output | 0.0023377 | 0.0052347 | 0.0081317 | 4.0 | 0.06 +Modify | 0.8194 | 0.82122 | 0.82303 | 0.2 | 9.20 +Other | | 0.06029 | | | 0.68 + +Nlocal: 51 ave 51 max 51 min +Histogram: 2 0 0 0 0 0 0 0 0 0 +Nghost: 3317 ave 3317 max 3317 min +Histogram: 2 0 0 0 0 0 0 0 0 0 +Neighs: 18320 ave 21941 max 14699 min +Histogram: 1 0 0 0 0 0 0 0 0 1 + +Total # of neighbors = 36640 +Ave neighs/atom = 359.21569 +Ave special neighs/atom = 2 +Neighbor list builds = 181 +Dangerous builds = 0 +Total wall time: 0:00:08 diff --git a/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.1 b/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.1 new file mode 100644 index 0000000000..a76489ebb3 --- /dev/null +++ b/examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.1 @@ -0,0 +1,437 @@ +LAMMPS (8 Feb 2023) +Processor partition = 1 + using 1 OpenMP thread(s) per MPI task +# Example for an alchemical transformation of two water molecules into a hydronium and hydroxyl ion +# WARNING: This input is intended for demonstrating the method only, +# the force field parameters are mostly made up and NOT suitable for production simulations. + +# set up different names for two partitions +variable name world twowater twoions + +units real +atom_style full +atom_modify map array +region box block -5 5 -5 5 -5 5 +boundary p p p +create_box 2 box bond/types 2 angle/types 2 extra/bond/per/atom 3 extra/angle/per/atom 3 extra/special/per/atom 3 +Created orthogonal box = (-5 -5 -5) to (5 5 5) + 1 by 1 by 2 MPI processor grid + +mass 1 15.9994 +mass 2 1.008 + +pair_style lj/cut/coul/cut 10.0 +pair_coeff 1 1 0.1553 3.166 +pair_coeff 1 2 0.0 1.0 +pair_coeff 2 2 0.0 1.0 + +bond_style harmonic +bond_coeff * 1000.0 1.0 + +angle_style harmonic +angle_coeff * 100.0 109.47 + +molecule water h2o.mol +Read molecule template water: + 1 molecules + 0 fragments + 3 atoms with max type 2 + 2 bonds with max type 1 + 1 angles with max type 1 + 0 dihedrals with max type 0 + 0 impropers with max type 0 + +# create the two molecules we want to transform ... +create_atoms 0 single -2.0 0.0 0.0 mol water 453624 +Created 3 atoms + using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5) + create_atoms CPU = 0.003 seconds +create_atoms 0 single 2.0 0.0 0.0 mol water 767353 +Created 3 atoms + using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5) + create_atoms CPU = 0.000 seconds + +# ... and put them in a group +group transform id 1:6 +6 atoms in group transform + +# now fill the rest of the box with more water +create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33 +Created 96 atoms + using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5) + create_atoms CPU = 0.001 seconds + +# change topology and settings for the two states +# we cannot simply create a different topology directly or +# load a different data file because the order and position +# of all atoms must be maintained across both replica + +# we first have to remove all topology data in the transform group +delete_bonds transform bond 1 +System init for delete_bonds ... +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 12 + ghost atom cutoff = 12 + binsize = 6, bins = 2 2 2 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair lj/cut/coul/cut, perpetual + attributes: half, newton on + pair build: half/bin/newton + stencil: half/bin/3d + bin: standard +Deleting bonds ... + 68 total bonds, 64 turned on, 4 turned off + 34 total angles, 34 turned on, 0 turned off + 0 total dihedrals, 0 turned on, 0 turned off + 0 total impropers, 0 turned on, 0 turned off +delete_bonds transform angle 1 remove +System init for delete_bonds ... +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +Deleting bonds ... + 64 total bonds, 64 turned on, 0 turned off + 32 total angles, 32 turned on, 0 turned off + 0 total dihedrals, 0 turned on, 0 turned off + 0 total impropers, 0 turned on, 0 turned off + +# then generate different topologies for the two partitions. select by name. +if "${name} == twowater" then "create_bonds single/bond 2 1 2" "create_bonds single/bond 2 1 3" "create_bonds single/bond 2 4 5" "create_bonds single/bond 2 4 6" "create_bonds single/angle 2 2 1 3" "create_bonds single/angle 2 5 4 6" else "create_bonds single/bond 2 1 2" "create_bonds single/bond 2 3 4" "create_bonds single/bond 2 4 5" "create_bonds single/bond 2 4 6" "create_bonds single/angle 2 3 4 5" "create_bonds single/angle 2 5 4 6" "create_bonds single/angle 2 3 4 6" "set atom 1 charge -1.1354" "set atom 2 charge 0.1354" "set atom 3 charge 0.56775" "set atom 4 charge -0.70305" "set atom 5*6 charge 0.56775" +create_bonds single/bond 2 1 2 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.002 seconds +create_bonds single/bond 2 3 4 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/bond 2 4 5 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 2 = max # of 1-2 neighbors + 1 = max # of 1-3 neighbors + 1 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/bond 2 4 6 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 3 = max # of 1-2 neighbors + 2 = max # of 1-3 neighbors + 2 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/angle 2 3 4 5 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 3 = max # of 1-2 neighbors + 2 = max # of 1-3 neighbors + 2 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/angle 2 5 4 6 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 3 = max # of 1-2 neighbors + 2 = max # of 1-3 neighbors + 2 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +create_bonds single/angle 2 3 4 6 +Finding 1-2 1-3 1-4 neighbors ... + special bond factors lj: 0 0 0 + special bond factors coul: 0 0 0 + 3 = max # of 1-2 neighbors + 2 = max # of 1-3 neighbors + 2 = max # of 1-4 neighbors + 7 = max # of special neighbors + special bonds CPU = 0.000 seconds +set atom 1 charge -1.1354 +Setting atom values ... + 1 settings made for charge +set atom 2 charge 0.1354 +Setting atom values ... + 1 settings made for charge +set atom 3 charge 0.56775 +Setting atom values ... + 1 settings made for charge +set atom 4 charge -0.70305 +Setting atom values ... + 1 settings made for charge +set atom 5*6 charge 0.56775 +Setting atom values ... + 2 settings made for charge + +velocity all create 300.0 5463576 +timestep 0.2 + +# define ramp variable to combine the two different partitions +if "${name} == twowater" then "variable ramp equal ramp(1.0,0.0)" else "variable ramp equal ramp(0.0,1.0)" +variable ramp equal ramp(0.0,1.0) + +# since the trajectory and forces are kept identical through fix alchemy, +# we can do fix npt simulations, but we must use the "mixed" pressure + +fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0 +fix transform all alchemy v_ramp +compute pressure all pressure/alchemy transform +fix_modify integrate press pressure + +# only need to output a dump file from one partition +# if "${name} == twowater" then # "dump 1 all atom 100 ${name}.lammpstrj" # "dump_modify 1 sort id" + +thermo_style custom step temp press etotal density pe ke f_transform f_transform[3] +thermo_modify colname f_transform lambda colname f_transform[3] EPot_mixed +thermo_modify press pressure + +thermo 100 +run 20000 +Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule +WARNING: Bond/angle/dihedral extent > half of periodic box length (src/domain.cpp:936) +Per MPI rank memory allocation (min/avg/max) = 7.535 | 7.535 | 7.535 Mbytes + Step Temp Press TotEng Density PotEng KinEng lambda EPot_mixed + 0 300 184631.17 10769.267 1.0171193 10678.948 90.318545 0 -892.51464 + 100 293.4335 -513.97927 21497.846 0.3029684 21409.505 88.341623 0.005 207.29243 + 200 320.81838 1513.3 3300.6279 0.64509407 3204.0418 96.586165 0.01 46.402434 + 300 261.53825 -3808.6942 -362.14861 0.88817871 -440.88779 78.739179 0.015 -614.76481 + 400 308.32444 1375.0327 -821.01681 0.97694485 -913.84152 92.824715 0.02 -897.24386 + 500 292.60623 -12843.384 -831.87625 1.0140814 -919.96881 88.092563 0.025 -1221.948 + 600 326.16112 91608.645 2941.9339 1.0263567 2843.7393 98.194659 0.03 2786.0743 + 700 317.22811 -60420.408 -1637.1259 0.9919143 -1732.6311 95.505272 0.035 -1616.7677 + 800 309.94673 -56793.931 -1501.6195 1.0024958 -1594.9326 93.313126 0.04 -1640.8606 + 900 272.80718 -5076.2282 -1452.719 1.0141078 -1534.8508 82.131824 0.045 -1357.1812 + 1000 346.52557 -16919.192 -1469.4601 1.0170322 -1573.7857 104.32562 0.05 -1514.1398 + 1100 300.62984 -51826.972 -2581.5818 1.0063592 -2672.09 90.508164 0.055 -2565.9668 + 1200 293.72767 45222.698 2045.9381 1.0272336 1957.5079 88.430186 0.06 1951.1683 + 1300 289.4971 94600.203 2042.1709 1.0298091 1955.0144 87.156524 0.065 2088.1011 + 1400 328.98666 35468.53 1805.9965 1.0198467 1706.9512 99.045321 0.07 1716.9268 + 1500 392.9407 41220.21 1663.5291 1.0194044 1545.2297 118.29944 0.075 1690.6416 + 1600 284.90093 -37995.141 -2488.5372 1.0068497 -2574.31 85.772791 0.08 -2529.5874 + 1700 267.67454 25275.255 1482.7666 1.0313556 1402.18 80.586582 0.085 1520.9695 + 1800 247.16495 -68222.994 -2869.7712 1.0018048 -2944.1832 74.41193 0.09 -2788.9545 + 1900 328.41849 -51979.665 -2578.9515 1.0154952 -2677.8258 98.874268 0.095 -2461.3422 + 2000 313.8715 21726.003 1310.0519 1.018831 1215.5571 94.494725 0.1 1363.0739 + 2100 296.07341 66130.11 1209.3035 1.0224882 1120.1671 89.136398 0.105 1272.1884 + 2200 282.67889 51746.464 2238.8823 1.0362683 2153.7785 85.103821 0.11 2276.1488 + 2300 281.04207 -71937.759 -2935.0422 1.0111968 -3019.6532 84.611035 0.115 -2892.2526 + 2400 321.00806 -51658.758 -2663.4071 1.0168161 -2760.0503 96.643269 0.12 -2622.9732 + 2500 288.95906 49645.753 2369.0667 1.0355648 2282.0722 86.99454 0.125 2275.0348 + 2600 303.0064 -48470.216 -2107.8769 1.000439 -2199.1005 91.223657 0.13 -2082.5895 + 2700 361.7083 54973.584 1906.9779 1.0223169 1798.0813 108.89656 0.135 1817.1972 + 2800 289.27827 33564.017 1687.0255 1.0257088 1599.9349 87.090642 0.14 1682.9162 + 2900 283.94365 -47190.274 -2449.1599 1.0059224 -2534.6445 85.484591 0.145 -2395.9575 + 3000 337.10996 61988.686 2157.6059 1.0258 2056.115 101.49094 0.15 2168.1085 + 3100 305.28249 -50123.109 -2137.6015 1.0143567 -2229.5104 91.9089 0.155 -2040.3453 + 3200 270.19471 -18114.265 -2175.7381 1.0132413 -2257.0834 81.345309 0.16 -2076.9595 + 3300 263.39996 -109301.89 -2594.4684 0.99956867 -2673.7681 79.299671 0.165 -2471.0269 + 3400 249.01722 -37585.949 -2655.0291 1.0080932 -2729.9987 74.969576 0.17 -2537.734 + 3500 267.88868 65528.865 2038.1308 1.0390872 1957.4797 80.651054 0.175 2061.6421 + 3600 351.46458 -67583.574 -3068.3126 1.0059621 -3174.1251 105.81256 0.18 -3004.5028 + 3700 332.82565 -61204.843 -2491.3087 1.0053485 -2591.5097 100.20109 0.185 -2415.4435 + 3800 298.0279 61123.085 1795.9521 1.020672 1706.2273 89.724819 0.19 1831.9734 + 3900 317.04715 43210.045 1685.6151 1.0287838 1590.1643 95.450792 0.195 1721.1961 + 4000 261.75174 -111920.72 -3478.9188 0.99457106 -3557.7223 78.803455 0.2 -3394.3528 + 4100 338.34245 89763.083 1919.3836 1.0318272 1817.5216 101.86199 0.205 1973.419 + 4200 288.52623 19989.792 1576.8543 1.0204119 1489.99 86.86423 0.21 1633.1439 + 4300 322.05593 -33791.403 -2303.9654 1.0155065 -2400.9242 96.958741 0.215 -2246.4353 + 4400 280.59054 -27839.599 -1776.7852 1.0123506 -1861.2603 84.475096 0.22 -1672.1425 + 4500 296.08336 75206.112 2655.8862 1.0329133 2566.7468 89.139393 0.225 2670.4908 + 4600 288.02094 -67445.101 -3254.8238 0.99579107 -3341.5359 86.712106 0.23 -3132.6352 + 4700 332.21641 -83964.647 -3046.3841 1.0053915 -3146.4017 100.01768 0.235 -2916.6652 + 4800 305.52557 -93111.682 -2635.4926 0.99824941 -2727.4747 91.982082 0.24 -2514.039 + 4900 272.82262 -53664.198 -3095.9932 0.99527907 -3178.1297 82.136472 0.245 -2947.2062 + 5000 288.96388 39846.089 1754.1018 1.0371205 1667.1058 86.99599 0.25 1833.3716 + 5100 293.25543 -71819.425 -3170.1841 0.99452016 -3258.4721 88.288011 0.255 -3069.7618 + 5200 296.16486 -64484.821 -2740.9832 1.0153125 -2830.1471 89.163931 0.26 -2589.6878 + 5300 342.07246 60255.985 1605.1961 1.0306374 1502.2111 102.98496 0.265 1608.5673 + 5400 295.53315 61421.189 1474.0782 1.0334762 1385.1045 88.973746 0.27 1492.0969 + 5500 285.1755 24845.113 981.19505 1.0335035 895.3396 85.855455 0.275 1014.0263 + 5600 276.79246 66244.072 1347.9794 1.0251134 1264.6478 83.331641 0.28 1390.9078 + 5700 288.43485 14629.129 1042.0147 1.02997 955.17799 86.836719 0.285 1074.3779 + 5800 276.64233 59262.97 680.78976 1.0332305 597.50331 83.286443 0.29 763.27415 + 5900 358.25384 -111203.23 -3652.2742 1.006331 -3760.1307 107.85655 0.295 -3561.6283 + 6000 327.82233 -62936.335 -2155.1414 1.0106476 -2253.8362 98.694786 0.3 -2062.107 + 6100 309.86553 -59110.381 -2958.2613 1.0076375 -3051.5499 93.28868 0.305 -2870.9566 + 6200 294.95752 -58011.845 -2023.2564 1.0161986 -2112.0569 88.800448 0.31 -1971.2586 + 6300 307.27252 59016.127 993.68575 1.0225343 901.17773 92.508022 0.315 1032.9616 + 6400 292.21611 -118857.83 -3289.8257 1.0031304 -3377.8009 87.975114 0.32 -3175.9068 + 6500 332.21108 97132.547 2078.6137 1.0349752 1978.5976 100.01607 0.325 2121.4674 + 6600 321.93985 -50017.836 -2763.8832 1.0076385 -2860.807 96.923796 0.33 -2673.6819 + 6700 299.07231 -98136.241 -3284.4174 1.0106261 -3374.4566 90.039252 0.335 -3194.0694 + 6800 285.33731 28487.407 232.87902 1.017186 146.97485 85.904169 0.34 294.02418 + 6900 296.41387 50054.533 1773.8315 1.0445934 1684.5926 89.238897 0.345 1851.1397 + 7000 318.34109 -71579.897 -3387.4573 1.007304 -3483.2977 95.840345 0.35 -3305.5715 + 7100 317.2192 12370.315 907.13585 1.0177238 811.63326 95.502589 0.355 961.68146 + 7200 289.91091 43944.333 972.20602 1.0302599 884.92492 87.281105 0.36 1039.3329 + 7300 324.81227 42859.173 994.92755 1.0300698 897.13897 97.788572 0.365 1083.5764 + 7400 305.72634 76516.72 2149.4427 1.0580673 2057.4002 92.042527 0.37 2218.3674 + 7500 283.77479 -87225.247 -3654.8517 0.9936609 -3740.2855 85.433754 0.375 -3581.9003 + 7600 318.72048 16084.995 1588.1487 1.0296163 1492.1941 95.954566 0.38 1653.5125 + 7700 307.36895 74007.024 1522.4279 1.0363553 1429.8909 92.537054 0.385 1572.8568 + 7800 297.44704 28693.429 1060.4277 1.0206402 970.87774 89.549946 0.39 1093.5903 + 7900 288.47173 -94644.76 -2723.8234 1.0095244 -2810.6712 86.847823 0.395 -2648.2974 + 8000 280.5001 63255.987 1380.8125 1.020587 1296.3647 84.44787 0.4 1418.7144 + 8100 320.23351 68383.584 2071.8202 1.0419027 1975.4101 96.410082 0.405 2111.2616 + 8200 354.86765 50012.316 1119.9483 1.0292074 1013.1112 106.8371 0.41 1149.0631 + 8300 278.8859 24200.145 486.00653 1.0232468 402.04464 83.961894 0.415 569.83331 + 8400 295.95754 -12010.401 800.25292 1.0200526 711.1514 89.101514 0.42 872.63254 + 8500 289.57525 -64114.236 -3307.2939 1.0062613 -3394.4739 87.18005 0.425 -3227.8592 + 8600 329.23192 15948.893 1155.4486 1.0271511 1056.3294 99.119161 0.43 1211.3868 + 8700 327.45998 39407.811 1425.7837 1.0343756 1327.198 98.585695 0.435 1480.9482 + 8800 296.99753 18017.036 587.91493 1.0178295 498.50031 89.414616 0.44 657.09357 + 8900 289.10016 79527.35 2286.7077 1.0462904 2199.6707 87.037018 0.445 2338.4742 + 9000 260.53308 -51430.717 -2982.8072 1.0040128 -3061.2437 78.436561 0.45 -2894.7913 + 9100 314.17867 40043.249 1796.3896 1.0361707 1701.8024 94.5872 0.455 1836.5399 + 9200 362.21236 -55153.163 -2572.3893 1.0102652 -2681.4376 109.04831 0.46 -2497.8765 + 9300 266.62049 -72193.959 -3333.755 1.0063615 -3414.0242 80.269248 0.465 -3225.2604 + 9400 259.59316 7674.9295 1105.0933 1.0326879 1026.9397 78.153588 0.47 1171.7119 + 9500 278.72962 -33127.835 -2875.2917 1.0163553 -2959.2065 83.914845 0.475 -2805.1711 + 9600 299.09039 24570.664 1594.6533 1.025317 1504.6086 90.044695 0.48 1631.8392 + 9700 336.65655 37767.059 1217.6584 1.0245842 1116.3039 101.35443 0.485 1218.4022 + 9800 279.68112 -55007.939 -2840.321 1.0129837 -2924.5223 84.201305 0.49 -2775.6909 + 9900 300.55694 14964.424 1101.8362 1.0176315 1011.35 90.486218 0.495 1127.7276 + 10000 276.77193 -75638.817 -3415.9268 1.0019055 -3499.2523 83.325459 0.5 -3341.663 + 10100 316.26926 61548.301 2199.0502 1.0427988 2103.8336 95.216598 0.505 2206.7159 + 10200 297.28171 56814.839 1828.6916 1.0270036 1739.1915 89.500171 0.51 1859.0176 + 10300 276.84774 83572.576 1662.6501 1.0322243 1579.3018 83.348283 0.515 1704.2233 + 10400 263.28687 -16549.994 903.97712 1.0177296 824.71149 79.265623 0.52 940.36339 + 10500 332.79252 47379.311 1277.3892 1.0223403 1177.1981 100.19112 0.525 1301.3762 + 10600 329.01105 39299.248 1039.7181 1.0230094 940.66548 99.052663 0.53 1077.5639 + 10700 339.21335 -4217.0313 724.12043 1.0235915 621.99624 102.12419 0.535 741.87272 + 10800 297.11192 -80013.424 -4028.1643 1.0009566 -4117.6134 89.449054 0.54 -3958.8283 + 10900 299.63096 47783.864 1609.2512 1.0435283 1519.0437 90.207441 0.545 1620.9328 + 11000 363.58243 67370.65 1557.8148 1.0450565 1448.354 109.46079 0.55 1543.803 + 11100 306.04664 50246.004 1020.8953 1.0341876 928.75637 92.138957 0.555 1063.5345 + 11200 255.68488 43750.967 1952.9201 1.0530607 1875.9431 76.976955 0.56 1986.9636 + 11300 271.71957 -40289.299 -3400.121 1.0168497 -3481.9254 81.804388 0.565 -3344.1558 + 11400 291.74007 -115249.24 -3603.3881 1.0067921 -3691.2199 87.831795 0.57 -3569.0256 + 11500 270.22077 28235.119 350.43952 1.0195886 269.08637 81.353155 0.575 373.78255 + 11600 297.08506 29790.701 793.45124 1.0317561 704.01027 89.440969 0.58 803.92835 + 11700 295.33119 -101140.58 -3395.5279 1.0107252 -3484.4409 88.912943 0.585 -3348.5652 + 11800 308.77287 -120312.75 -3874.1096 0.99681611 -3967.0693 92.95972 0.59 -3835.4663 + 11900 293.86337 82287.969 1623.3312 1.0331446 1534.8602 88.47104 0.595 1625.4846 + 12000 279.9866 27919.689 2021.9956 1.0384909 1937.7024 84.293273 0.6 2028.6311 + 12100 292.19439 47035.018 962.98803 1.0210948 875.01946 87.968573 0.605 976.68113 + 12200 317.78325 46886.483 1744.1855 1.0310567 1648.5131 95.672401 0.61 1763.6359 + 12300 357.44038 51492.35 1650.3113 1.0414835 1542.6996 107.61165 0.615 1638.1905 + 12400 288.16042 -50528.478 -3145.341 1.0130798 -3232.0951 86.7541 0.62 -3115.1383 + 12500 300.37356 -114034.29 -3383.5763 1.0130272 -3474.0074 90.431009 0.625 -3364.511 + 12600 265.4402 73960.28 1102.3119 1.0343798 1022.398 79.913907 0.63 1094.544 + 12700 287.07994 22722.337 1499.2541 1.0268893 1412.8252 86.428809 0.635 1490.443 + 12800 315.39745 46072.191 663.05431 1.0205906 568.10018 94.954129 0.64 662.83753 + 12900 307.48549 -120472.73 -3484.9054 1.005689 -3577.4775 92.57214 0.645 -3469.7834 + 13000 291.00577 -59154.878 -3000.6461 1.012821 -3088.2569 87.610727 0.65 -2966.2448 + 13100 318.37718 86703.775 1825.1817 1.036209 1729.3305 95.851211 0.655 1819.7078 + 13200 280.41509 -107033.79 -2964.9455 1.0135499 -3049.3678 84.422277 0.66 -2953.8355 + 13300 275.01919 112086.85 2399.5628 1.0437945 2316.765 82.797775 0.665 2394.6707 + 13400 299.46547 -111959.62 -3814.3709 0.98732074 -3904.5286 90.157617 0.67 -3792.6798 + 13500 275.65631 17625.49 1528.0027 1.0255473 1445.0131 82.989588 0.675 1521.2745 + 13600 302.63525 122466.63 2631.7819 1.0438137 2540.67 91.111917 0.68 2611.7188 + 13700 302.21323 -94399.663 -3344.9423 1.0068914 -3435.9271 90.984863 0.685 -3328.4539 + 13800 305.70483 30039.825 1506.1765 1.0225478 1414.1404 92.036052 0.69 1502.8476 + 13900 273.15564 -27822.597 -2713.3946 1.0131002 -2795.6313 82.236732 0.695 -2683.7804 + 14000 257.14635 -1821.0012 1045.462 1.0172653 968.04508 77.416947 0.7 1067.8821 + 14100 323.74174 -76139.363 -3042.0486 1.0139806 -3139.5149 97.466276 0.705 -3052.598 + 14200 274.99913 54440.943 796.17718 1.0186356 713.38545 82.791736 0.71 767.64274 + 14300 291.8774 1281.7286 1079.6881 1.0214356 991.81498 87.873139 0.715 1042.8585 + 14400 308.97245 -47961.551 -3005.1106 1.0155221 -3098.1304 93.019806 0.72 -3007.3204 + 14500 311.05309 86026.739 1851.7257 1.0311818 1758.0795 93.646208 0.725 1830.5463 + 14600 269.27432 -12972.747 726.82371 1.0172239 645.7555 81.068214 0.73 711.81484 + 14700 262.20739 -77722.459 -3545.7104 1.0029421 -3624.651 78.940634 0.735 -3537.1929 + 14800 300.33592 -107332.8 -3118.6094 1.0126785 -3209.0291 90.419677 0.74 -3125.8439 + 14900 307.0173 -73204.94 -3016.2921 1.016875 -3108.7233 92.431185 0.745 -3016.0143 + 15000 334.08711 64440.709 1360.8856 1.0369481 1260.3047 100.58087 0.75 1323.9179 + 15100 293.46703 35179.395 1804.4648 1.0331703 1716.1131 88.351717 0.755 1768.9681 + 15200 307.80546 54006.291 852.72781 1.0246134 760.05934 92.668471 0.76 822.22754 + 15300 336.59792 33832.907 988.39905 1.0242407 887.06227 101.33678 0.765 956.25285 + 15400 285.46526 -99629.94 -2736.0944 1.0041261 -2822.0371 85.942688 0.77 -2730.1887 + 15500 269.15938 122546.45 2708.063 1.0318941 2627.0294 81.033612 0.775 2695.4331 + 15600 294.50818 -104733.19 -2271.5729 1.0017241 -2360.238 88.665166 0.78 -2281.5773 + 15700 321.05125 -40438.501 -2060.144 1.0125007 -2156.8002 96.656271 0.785 -2084.0893 + 15800 303.22185 -42224.408 -2486.7037 1.0090867 -2577.9922 91.288521 0.79 -2506.3355 + 15900 281.65025 78652.063 2532.2753 1.0460521 2447.4812 84.794136 0.795 2500.7128 + 16000 286.5056 63174.412 1479.3631 1.0185277 1393.1072 86.255897 0.8 1442.2427 + 16100 271.74574 -75051.584 -2429.2879 1.0091727 -2511.1002 81.812265 0.805 -2447.1967 + 16200 266.40007 -64120.09 -2497.2609 1.0064388 -2577.4638 80.202888 0.81 -2525.4243 + 16300 330.67137 46204.124 1582.1688 1.0271037 1482.6162 99.552523 0.815 1529.8186 + 16400 310.94573 18510.676 1235.5974 1.0188695 1141.9835 93.613887 0.82 1195.3041 + 16500 270.29191 79005.409 1711.3568 1.0354088 1629.9822 81.374574 0.825 1678.7342 + 16600 334.24108 17143.221 1446.3175 1.0230561 1345.6903 100.62723 0.83 1379.3215 + 16700 292.28721 -78618.031 -3045.0258 0.99940789 -3133.0223 87.996518 0.835 -3087.379 + 16800 346.09157 -69398.937 -2956.6793 1.0153026 -3060.8742 104.19496 0.84 -3011.3476 + 16900 284.37042 39314.526 2102.7814 1.0343186 2017.1683 85.613074 0.845 2045.8316 + 17000 281.69389 86032.163 1577.5995 1.0237583 1492.7922 84.807275 0.85 1528.8391 + 17100 297.2629 55128.393 2431.5835 1.0419182 2342.089 89.494508 0.855 2369.0347 + 17200 328.81472 47063.86 1584.1492 1.0364627 1485.1556 98.993555 0.86 1526.7115 + 17300 313.41363 -46685.063 -2635.2463 1.008117 -2729.6032 94.356878 0.865 -2685.8034 + 17400 288.17114 68721.523 2548.303 1.0429105 2461.5456 86.757326 0.87 2490.5845 + 17500 260.01702 -23324.539 -1776.7039 1.0100415 -1854.9851 78.281195 0.875 -1812.4813 + 17600 297.6633 72795.463 2396.9788 1.0362827 2307.3638 89.615054 0.88 2336.0152 + 17700 301.02937 93833.118 2424.9649 1.0276466 2334.3364 90.628448 0.885 2364.8751 + 17800 238.92837 -78317.212 -2548.5755 1.0084077 -2620.5077 71.932208 0.89 -2579.4429 + 17900 299.27724 -65278.891 -2933.564 1.0034773 -3023.6649 90.10095 0.895 -2987.498 + 18000 315.15753 -72437.481 -2654.8693 1.0136139 -2749.7512 94.881898 0.9 -2718.125 + 18100 268.39411 -42194.209 -2606.045 1.0108104 -2686.8482 80.803219 0.905 -2656.3853 + 18200 284.91941 40905.162 1786.1505 1.0266185 1700.3722 85.778354 0.91 1723.6803 + 18300 329.41629 -56526.814 -2472.7835 1.0137784 -2571.9582 99.174667 0.915 -2540.6829 + 18400 293.16732 -22085.551 -1991.6537 1.0163557 -2079.9152 88.261484 0.92 -2052.7864 + 18500 276.28141 29367.359 2034.0833 1.0270087 1950.9056 83.177784 0.925 1971.9405 + 18600 331.18031 94669.124 2160.1352 1.0225403 2060.4295 99.705745 0.93 2081.8675 + 18700 265.4491 49872.536 2324.2391 1.0255308 2244.3225 79.916587 0.935 2260.4261 + 18800 314.51042 60758.629 2298.6407 1.0300818 2203.9536 94.687078 0.94 2219.6547 + 18900 317.70309 66345.835 1490.6939 1.0248229 1395.0457 95.648268 0.945 1410.3087 + 19000 291.39831 -96336.931 -2732.3994 1.0100498 -2820.1283 87.728905 0.95 -2805.9396 + 19100 303.37353 -21421.495 -2360.5644 0.99771284 -2451.8986 91.334186 0.955 -2436.2173 + 19200 293.92948 36383.881 1708.7388 1.0251299 1620.2479 88.490944 0.96 1633.579 + 19300 287.24409 -46691.027 -2736.6663 1.0112067 -2823.1445 86.478228 0.965 -2811.3506 + 19400 333.42205 21152.63 1209.1556 1.0189912 1108.7749 100.38065 0.97 1118.2463 + 19500 309.17989 60445.761 2014.5979 1.0338543 1921.5157 93.082259 0.975 1928.6181 + 19600 255.92246 -32795.547 -2145.3103 1.0157044 -2222.3587 77.04848 0.98 -2213.9484 + 19700 342.62332 33538.746 2181.2002 1.0225905 2078.0494 103.1508 0.985 2081.9185 + 19800 295.68946 81474.342 1996.4567 1.0240743 1907.4359 89.020804 0.99 1910.5545 + 19900 306.4947 43488.052 2048.2997 1.0222707 1956.0258 92.27385 0.995 1957.4368 + 20000 313.31679 -25133.284 -1555.4274 1.0163289 -1649.7551 94.327722 1 -1649.7551 +Loop time of 8.92592 on 2 procs for 20000 steps with 102 atoms + +Performance: 38.719 ns/day, 0.620 hours/ns, 2240.664 timesteps/s, 228.548 katom-step/s +96.2% CPU use with 2 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.7753 | 6.0676 | 7.3599 | 52.5 | 67.98 +Bond | 0.036352 | 0.036909 | 0.037466 | 0.3 | 0.41 +Neigh | 0.085317 | 0.085383 | 0.085449 | 0.0 | 0.96 +Comm | 0.72794 | 2.022 | 3.3161 | 91.0 | 22.65 +Output | 0.0023255 | 0.005156 | 0.0079865 | 3.9 | 0.06 +Modify | 0.64567 | 0.64689 | 0.64812 | 0.2 | 7.25 +Other | | 0.06196 | | | 0.69 + +Nlocal: 51 ave 51 max 51 min +Histogram: 2 0 0 0 0 0 0 0 0 0 +Nghost: 3317 ave 3317 max 3317 min +Histogram: 2 0 0 0 0 0 0 0 0 0 +Neighs: 18319.5 ave 21940 max 14699 min +Histogram: 1 0 0 0 0 0 0 0 0 1 + +Total # of neighbors = 36639 +Ave neighs/atom = 359.20588 +Ave special neighs/atom = 2.0196078 +Neighbor list builds = 181 +Dangerous builds = 0 +Total wall time: 0:00:08 diff --git a/src/.gitignore b/src/.gitignore index 8f0f577a45..0d73a7911d 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1465,6 +1465,8 @@ /zstd_file_writer.cpp /zstd_file_writer.h +/compute_pressure_alchemy.cpp +/compute_pressure_alchemy.h /atom_vec_smd.cpp /atom_vec_smd.h /compute_saed.cpp @@ -1531,8 +1533,8 @@ /fix_langevin_drude.h /fix_mol_swap.cpp /fix_mol_swap.h -/fix_pimd.cpp -/fix_pimd.h +/fix_alchemy.cpp +/fix_alchemy.h /fix_pimd_nvt.cpp /fix_pimd_nvt.h /fix_qbmsst.cpp diff --git a/src/REPLICA/compute_pressure_alchemy.cpp b/src/REPLICA/compute_pressure_alchemy.cpp new file mode 100644 index 0000000000..d19ea1b893 --- /dev/null +++ b/src/REPLICA/compute_pressure_alchemy.cpp @@ -0,0 +1,102 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "compute_pressure_alchemy.h" + +#include "domain.h" +#include "error.h" +#include "fix.h" +#include "modify.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputePressureAlchemy::ComputePressureAlchemy(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR, "Illegal compute pressure/alchemy command"); + if (igroup) error->all(FLERR, "Compute pressure/alchemy must use group all"); + + scalar_flag = vector_flag = 1; + size_vector = 6; + extscalar = 0; + extvector = 0; + pressflag = 1; + timeflag = 1; + + id_fix = arg[3]; + if (!modify->get_fix_by_id(id_fix)) + error->all(FLERR, "Could not find compute pressure/alchemy fix ID {} for fix alchemy", id_fix); + + vector = new double[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputePressureAlchemy::~ComputePressureAlchemy() +{ + delete[] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputePressureAlchemy::init() +{ + + fix = modify->get_fix_by_id(id_fix); + if (!fix) + error->all(FLERR, "Could not find compute pressure/alchemy fix ID {} for fix alchemy", id_fix); + + int dim = 0; + void *ptr = fix->extract("pressure", dim); + if (!ptr || (dim != 1)) error->all(FLERR, "Could not extract pressure from fix alchemy"); +} + +/* ---------------------------------------------------------------------- + compute total pressure from tensor, averaged over Pxx, Pyy, Pzz +------------------------------------------------------------------------- */ + +double ComputePressureAlchemy::compute_scalar() +{ + invoked_scalar = update->ntimestep; + if (update->vflag_global != invoked_scalar) + error->all(FLERR, "Virial was not tallied on needed timestep"); + + compute_vector(); + + if (domain->dimension == 3) { + scalar = (vector[0] + vector[1] + vector[2]) / 3.0; + } else { + scalar = (vector[0] + vector[1]) / 2.0; + } + return scalar; +} + +/* ---------------------------------------------------------------------- + extract compute combined system pressure tensor from alchemy fix +------------------------------------------------------------------------- */ + +void ComputePressureAlchemy::compute_vector() +{ + invoked_vector = update->ntimestep; + if (update->vflag_global != invoked_vector) + error->all(FLERR, "Virial was not tallied on needed timestep"); + + int dim = 0; + double *pressure = (double *) fix->extract("pressure", dim); + if (!pressure || (dim != 1)) error->all(FLERR, "Could not extract pressure from fix alchemy"); + + for (int i = 0; i < 6; i++) vector[i] = pressure[i]; +} diff --git a/src/REPLICA/compute_pressure_alchemy.h b/src/REPLICA/compute_pressure_alchemy.h new file mode 100644 index 0000000000..0ce0039587 --- /dev/null +++ b/src/REPLICA/compute_pressure_alchemy.h @@ -0,0 +1,41 @@ +/* -*- 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 COMPUTE_CLASS +// clang-format off +ComputeStyle(pressure/alchemy,ComputePressureAlchemy); +// clang-format on +#else + +#ifndef LMP_COMPUTE_PRESSURE_ALCHEMY_H +#define LMP_COMPUTE_PRESSURE_ALCHEMY_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputePressureAlchemy : public Compute { + public: + ComputePressureAlchemy(class LAMMPS *, int, char **); + ~ComputePressureAlchemy() override; + void init() override; + double compute_scalar() override; + void compute_vector() override; + + protected: + class Fix *fix; + std::string id_fix; +}; +} // namespace LAMMPS_NS +#endif +#endif diff --git a/src/REPLICA/fix_alchemy.cpp b/src/REPLICA/fix_alchemy.cpp new file mode 100644 index 0000000000..2f99f667b6 --- /dev/null +++ b/src/REPLICA/fix_alchemy.cpp @@ -0,0 +1,265 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "fix_alchemy.h" + +#include "atom.h" +#include "comm.h" +#include "compute.h" +#include "domain.h" +#include "error.h" +#include "input.h" +#include "memory.h" +#include "modify.h" +#include "respa.h" +#include "universe.h" +#include "update.h" +#include "variable.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), commbuf(nullptr) +{ + if (narg != 4) error->all(FLERR, "Incorrect number of arguments for fix alchemy"); + if (universe->nworlds != 2) error->all(FLERR, "Must use exactly two partitions"); + if (utils::strmatch(arg[3], "^v_")) + id_lambda = arg[3] + 2; + else + error->all(FLERR, "Must use variable as lambda argument to fix alchemy"); + + lambda = epot[0] = epot[1] = epot[2] = 0.0; + progress = 0; + for (int i = 0; i < 6; ++i) pressure[i] = 0.0; + + no_change_box = 1; + time_depend = 1; + scalar_flag = 1; + extscalar = 0; + vector_flag = 1; + size_vector = 3; + extvector = 1; + ilevel_respa = 0; + nmax = 6; + sync_box = 0; + + // set up rank-to-rank communicator for inter-partition communication + + int color = comm->me; + int key = universe->iworld; + MPI_Comm_split(universe->uworld, color, key, &samerank); + + // check that we have the same domain decomposition on all ranks + + int my_nlocal[2] = {0, 0}; + int all_nlocal[2] = {0, 0}; + my_nlocal[universe->iworld] = atom->nlocal; + MPI_Allreduce(my_nlocal, all_nlocal, 2, MPI_INT, MPI_SUM, samerank); + int fail = (all_nlocal[0] == all_nlocal[1]) ? 0 : 1; + int allfail = 0; + MPI_Allreduce(&fail, &allfail, 1, MPI_INT, MPI_MAX, universe->uworld); + if (allfail) + error->all(FLERR, "Number of atoms and domain decomposition must match for both partitions"); + + id_pe = std::string(id) + "_pe"; + pe = modify->add_compute(id_pe + " all pe"); + pe->addstep(update->ntimestep); + id_temp = std::string(id) + "_temp"; + temp = modify->add_compute(id_temp + " all temp"); + temp->addstep(update->ntimestep); + id_press = std::string(id) + "_press"; + press = modify->add_compute(id_press + " all pressure " + id_temp); + press->addstep(update->ntimestep); +} + +/* ---------------------------------------------------------------------- */ + +FixAlchemy::~FixAlchemy() +{ + MPI_Comm_free(&samerank); + modify->delete_compute(id_pe); + modify->delete_compute(id_temp); + modify->delete_compute(id_press); + memory->destroy(commbuf); +} + +/* ---------------------------------------------------------------------- */ + +int FixAlchemy::setmask() +{ + int mask = 0; + mask |= POST_INTEGRATE; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +static void synchronize_box(Domain *domain, MPI_Comm samerank) +{ + MPI_Bcast(&domain->boxlo[0], 3, MPI_DOUBLE, 0, samerank); + MPI_Bcast(&domain->boxhi[0], 3, MPI_DOUBLE, 0, samerank); + MPI_Bcast(&domain->yz, 1, MPI_DOUBLE, 0, samerank); + MPI_Bcast(&domain->xz, 1, MPI_DOUBLE, 0, samerank); + MPI_Bcast(&domain->xy, 1, MPI_DOUBLE, 0, samerank); + domain->set_global_box(); + domain->set_local_box(); +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::init() +{ + int onenmax = MAX(nmax, 3 * atom->nmax); + MPI_Allreduce(&onenmax, &nmax, 1, MPI_INT, MPI_MAX, universe->uworld); + memory->destroy(commbuf); + memory->create(commbuf, sizeof(double) * nmax, "alchemy:nmax"); + + if (modify->get_fix_by_style("^balance").size() > 0) + error->all(FLERR, "Fix alchemy is not compatible with load balancing"); + + if (modify->get_fix_by_style("^alchemy").size() > 1) + error->all(FLERR, "There may only one fix alchemy at a time"); + + ivar = input->variable->find(id_lambda.c_str()); + if (ivar < 0) + error->universe_one(FLERR, fmt::format("Variable {} for fix alchemy does not exist", id_lambda)); + if (!input->variable->equalstyle(ivar)) + error->universe_one(FLERR, + fmt::format("Variable {} for fix alchemy is invalid style", id_lambda)); + lambda = input->variable->compute_equal(ivar); + + // synchronize box dimensions, determine if resync during run will be needed. + + synchronize_box(domain, samerank); + + sync_box = 0; + for (auto ifix : modify->get_fix_list()) + if (ifix->box_change) sync_box = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::setup(int vflag) +{ + if (utils::strmatch(update->integrate_style, "^respa")) { + auto respa = dynamic_cast(update->integrate); + respa->copy_flevel_f(ilevel_respa); + post_force_respa(vflag, ilevel_respa, 0); + respa->copy_f_flevel(ilevel_respa); + } else { + post_force(vflag); + } + + if (universe->me == 0) { + double delta = update->ntimestep - update->beginstep; + if ((delta != 0.0) && (update->beginstep != update->endstep)) + delta /= update->endstep - update->beginstep; + progress = static_cast(delta*100.0); + auto msg = fmt::format("Starting alchemical transformation at {:>3d}%\n", progress); + if (universe->uscreen) fmt::print(universe->uscreen, msg); + if (universe->ulogfile) fmt::print(universe->ulogfile, msg); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::post_integrate() +{ + // synchronize atom positions + + const int nall = atom->nlocal + atom->nghost; + MPI_Bcast(&atom->x[0][0], 3 * nall, MPI_DOUBLE, 0, samerank); + + // synchronize box dimensions, if needed + + if (sync_box) synchronize_box(domain, samerank); +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::post_force(int /*vflag*/) +{ + if (3 * atom->nmax > nmax) { + nmax = 3 * atom->nmax; + memory->grow(commbuf, sizeof(double) * atom->nmax, "alchemy:commbuf"); + } + + const int nall = 3 * atom->nlocal; + double *f = &atom->f[0][0]; + lambda = input->variable->compute_equal(ivar); + + for (int i = 0; i < nall; ++i) commbuf[i] = f[i] * lambda; + MPI_Allreduce(commbuf, f, nall, MPI_DOUBLE, MPI_SUM, samerank); + + // sum up potential energy + + const double scalefac = 1.0 / comm->nprocs; + commbuf[0] = commbuf[1] = commbuf[2] = 0.0; + commbuf[universe->iworld] = scalefac * pe->compute_scalar(); + commbuf[2] = lambda * scalefac * pe->compute_scalar(); + MPI_Allreduce(commbuf, epot, 3, MPI_DOUBLE, MPI_SUM, universe->uworld); + pe->addstep(update->ntimestep + 1); + + // sum up pressure + + press->compute_vector(); + for (int i = 0; i < 6; ++i) commbuf[i] = lambda * scalefac * press->vector[i]; + MPI_Allreduce(commbuf, pressure, 6, MPI_DOUBLE, MPI_SUM, universe->uworld); + press->addstep(update->ntimestep + 1); + + // print progress info + + if (universe->me == 0) { + double delta = update->ntimestep - update->beginstep; + if ((delta != 0.0) && (update->beginstep != update->endstep)) + delta /= update->endstep - update->beginstep; + int status = static_cast(delta*100.0); + if ((status / 10) > (progress / 10)) { + progress = status; + auto msg = fmt::format(" Alchemical transformation progress: {:>3d}%\n", progress); + if (universe->uscreen) fmt::print(universe->uscreen, msg); + if (universe->ulogfile) fmt::print(universe->ulogfile, msg); + } + } +} + +/* ---------------------------------------------------------------------- */ + +double FixAlchemy::compute_scalar() +{ + return lambda; +} + +/* ---------------------------------------------------------------------- */ + +double FixAlchemy::compute_vector(int n) +{ + return epot[n]; +} + +/* ---------------------------------------------------------------------- */ + +void *FixAlchemy::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str, "lambda") == 0) { return λ } + if (strcmp(str, "pe") == 0) { return &epot[2]; } + dim = 1; + if (strcmp(str, "pressure") == 0) { return pressure; } + return nullptr; +} diff --git a/src/REPLICA/fix_alchemy.h b/src/REPLICA/fix_alchemy.h new file mode 100644 index 0000000000..677a29b643 --- /dev/null +++ b/src/REPLICA/fix_alchemy.h @@ -0,0 +1,58 @@ +/* -*- 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 FIX_CLASS +// clang-format off +FixStyle(alchemy,FixAlchemy); +// clang-format on +#else + +#ifndef LMP_FIX_ALCHEMY_H +#define LMP_FIX_ALCHEMY_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAlchemy : public Fix { + public: + FixAlchemy(class LAMMPS *, int, char **); + ~FixAlchemy() override; + + int setmask() override; + void init() override; + void setup(int) override; + void post_integrate() override; + void post_force(int) override; + double compute_scalar() override; + double compute_vector(int) override; + void *extract(const char *, int &) override; + + protected: + MPI_Comm samerank; + double *commbuf; + class Compute *pe, *temp, *press; + std::string id_pe, id_temp, id_press, id_lambda; + double lambda; // scaling prefactor for combining the partitions + double epot[3]; // last (unscaled) potential energy from each replica and combined energy + double pressure[6]; // joined pressure + int progress; // for progress indicator + int sync_box; // 1 of box dimensions need to be synchronized + int ilevel_respa; + int nmax; + int ivar; +}; +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/create_bonds.cpp b/src/create_bonds.cpp index c6d3f1838e..5355ce2b7c 100644 --- a/src/create_bonds.cpp +++ b/src/create_bonds.cpp @@ -319,8 +319,6 @@ void CreateBonds::many() void CreateBonds::single_bond() { - int m; - // check that 2 atoms exist const int nlocal = atom->nlocal; @@ -341,7 +339,8 @@ void CreateBonds::single_bond() int **bond_type = atom->bond_type; tagint **bond_atom = atom->bond_atom; - if ((m = idx1) >= 0) { + int m = idx1; + if ((m >= 0) && (m < nlocal)) { if (num_bond[m] == atom->bond_per_atom) error->one(FLERR, "New bond exceeded bonds per atom in create_bonds"); bond_type[m][num_bond[m]] = btype; @@ -352,7 +351,8 @@ void CreateBonds::single_bond() if (force->newton_bond) return; - if ((m = idx2) >= 0) { + m = idx2; + if ((m >= 0) && (m < nlocal)) { if (num_bond[m] == atom->bond_per_atom) error->one(FLERR, "New bond exceeded bonds per atom in create_bonds"); bond_type[m][num_bond[m]] = btype; @@ -365,8 +365,6 @@ void CreateBonds::single_bond() void CreateBonds::single_angle() { - int m; - // check that 3 atoms exist const int nlocal = atom->nlocal; @@ -391,7 +389,8 @@ void CreateBonds::single_angle() tagint **angle_atom2 = atom->angle_atom2; tagint **angle_atom3 = atom->angle_atom3; - if ((m = idx2) >= 0) { + int m = idx2; + if ((m >= 0) && (m < nlocal)) { if (num_angle[m] == atom->angle_per_atom) error->one(FLERR, "New angle exceeded angles per atom in create_bonds"); angle_type[m][num_angle[m]] = atype; @@ -404,7 +403,8 @@ void CreateBonds::single_angle() if (force->newton_bond) return; - if ((m = idx1) >= 0) { + m = idx1; + if ((m >= 0) && (m < nlocal)) { if (num_angle[m] == atom->angle_per_atom) error->one(FLERR, "New angle exceeded angles per atom in create_bonds"); angle_type[m][num_angle[m]] = atype; @@ -414,7 +414,8 @@ void CreateBonds::single_angle() num_angle[m]++; } - if ((m = idx3) >= 0) { + m = idx3; + if ((m >= 0) && (m < nlocal)) { if (num_angle[m] == atom->angle_per_atom) error->one(FLERR, "New angle exceeded angles per atom in create_bonds"); angle_type[m][num_angle[m]] = atype; @@ -429,8 +430,6 @@ void CreateBonds::single_angle() void CreateBonds::single_dihedral() { - int m; - // check that 4 atoms exist const int nlocal = atom->nlocal; @@ -458,7 +457,8 @@ void CreateBonds::single_dihedral() tagint **dihedral_atom3 = atom->dihedral_atom3; tagint **dihedral_atom4 = atom->dihedral_atom4; - if ((m = idx2) >= 0) { + int m = idx2; + if ((m >= 0) && (m < nlocal)) { if (num_dihedral[m] == atom->dihedral_per_atom) error->one(FLERR, "New dihedral exceeded dihedrals per atom in create_bonds"); dihedral_type[m][num_dihedral[m]] = dtype; @@ -472,7 +472,8 @@ void CreateBonds::single_dihedral() if (force->newton_bond) return; - if ((m = idx1) >= 0) { + m = idx1; + if ((m >= 0) && (m < nlocal)) { if (num_dihedral[m] == atom->dihedral_per_atom) error->one(FLERR, "New dihedral exceeded dihedrals per atom in create_bonds"); dihedral_type[m][num_dihedral[m]] = dtype; @@ -483,7 +484,8 @@ void CreateBonds::single_dihedral() num_dihedral[m]++; } - if ((m = idx3) >= 0) { + m = idx3; + if ((m >= 0) && (m < nlocal)) { if (num_dihedral[m] == atom->dihedral_per_atom) error->one(FLERR, "New dihedral exceeded dihedrals per atom in create_bonds"); dihedral_type[m][num_dihedral[m]] = dtype; @@ -494,7 +496,8 @@ void CreateBonds::single_dihedral() num_dihedral[m]++; } - if ((m = idx4) >= 0) { + m = idx4; + if ((m >= 0) && (m < nlocal)) { if (num_dihedral[m] == atom->dihedral_per_atom) error->one(FLERR, "New dihedral exceeded dihedrals per atom in create_bonds"); dihedral_type[m][num_dihedral[m]] = dtype; @@ -510,8 +513,6 @@ void CreateBonds::single_dihedral() void CreateBonds::single_improper() { - int m; - // check that 4 atoms exist const int nlocal = atom->nlocal; @@ -539,7 +540,8 @@ void CreateBonds::single_improper() tagint **improper_atom3 = atom->improper_atom3; tagint **improper_atom4 = atom->improper_atom4; - if ((m = idx2) >= 0) { + int m = idx2; + if ((m >= 0) && (m < nlocal)) { if (num_improper[m] == atom->improper_per_atom) error->one(FLERR, "New improper exceeded impropers per atom in create_bonds"); improper_type[m][num_improper[m]] = dtype; @@ -553,7 +555,8 @@ void CreateBonds::single_improper() if (force->newton_bond) return; - if ((m = idx1) >= 0) { + m = idx1; + if ((m >= 0) && (m < nlocal)) { if (num_improper[m] == atom->improper_per_atom) error->one(FLERR, "New improper exceeded impropers per atom in create_bonds"); improper_type[m][num_improper[m]] = dtype; @@ -564,7 +567,8 @@ void CreateBonds::single_improper() num_improper[m]++; } - if ((m = idx3) >= 0) { + m = idx3; + if ((m >= 0) && (m < nlocal)) { if (num_improper[m] == atom->improper_per_atom) error->one(FLERR, "New improper exceeded impropers per atom in create_bonds"); improper_type[m][num_improper[m]] = dtype; @@ -575,7 +579,8 @@ void CreateBonds::single_improper() num_improper[m]++; } - if ((m = idx4) >= 0) { + m = idx4; + if ((m >= 0) && (m < nlocal)) { if (num_improper[m] == atom->improper_per_atom) error->one(FLERR, "New improper exceeded impropers per atom in create_bonds"); improper_type[m][num_improper[m]] = dtype;