From 648cd3f0c6a7a54d3eb99173dab1e707d978674e Mon Sep 17 00:00:00 2001 From: talinke Date: Tue, 22 Apr 2025 17:43:21 -0700 Subject: [PATCH] Temperature bias capability --- examples/gjf/README | 32 +++--- examples/gjf/in.gjf.vfull | 4 +- examples/gjf/in.gjf.vhalf | 2 +- examples/gjf/log.2Apr25.gjf.vfull.g++.1 | 50 ++++---- examples/gjf/log.2Apr25.gjf.vfull.g++.4 | 50 ++++---- examples/gjf/log.2Apr25.gjf.vhalf.g++.1 | 59 ++++++---- examples/gjf/log.2Apr25.gjf.vhalf.g++.4 | 61 ++++++---- src/EXTRA-FIX/fix_langevin_gjf.cpp | 145 +++++++++++++++++------- src/EXTRA-FIX/fix_langevin_gjf.h | 2 +- 9 files changed, 242 insertions(+), 163 deletions(-) diff --git a/examples/gjf/README b/examples/gjf/README index b98e4134ab..34b60292bb 100644 --- a/examples/gjf/README +++ b/examples/gjf/README @@ -15,23 +15,23 @@ with the correct Boltzmann statistics. A comparison using averaged properties from this example's input file is shown below. 'X' denotes a failed simulation. The theoretical value for KE is 1.1168 eV. -KINETIC ENERGY (eV) -| Δt || 0.01 | 0.05 | 0.10 | 0.11 | 0.12 | 0.13 | 0.14 | -|==============||========|========|========|========|========|========|========| -| gjf half || 1.117 | 1.116 | 1.119 | 1.119 | 1.123 | 1.136 | 1.170 | -| gjf full || 1.116 | 1.071 | 0.938 | 0.898 | 0.858 | 0.817 | 0.780 | -| langevin || 1.110 | 1.113 | 1.121 | 1.129 | 1.157 | X | X | -| nvt || 1.083 | 1.109 | 1.112 | 1.113 | 1.114 | X | X | -|--------------||--------|--------|--------|--------|--------|--------|--------| - POTENTIAL ENERGY (eV) -| Δt || 0.01 | 0.05 | 0.10 | 0.11 | 0.12 | 0.13 | 0.14 | -|==============||========|========|========|========|========|========|========| -| gjf half || -55.11 | -55.11 | -55.11 | -55.11 | -55.11 | -55.10 | -55.07 | -| gjf full || -55.11 | -55.11 | -55.11 | -55.11 | -55.11 | -55.10 | -55.07 | -| langevin || -55.11 | -55.07 | -54.87 | -54.79 | -54.65 | X | X | -| nvt || -55.14 | -55.07 | -54.90 | -54.84 | -54.76 | X | X | -|--------------||--------|--------|--------|--------|--------|--------|--------| +| Δt || 0.01 | 0.05 | 0.10 | 0.11 | 0.12 | 0.13 | 0.14 | +|===================||========|========|========|========|========|========|========| +| gjf half || -55.11 | -55.11 | -55.11 | -55.11 | -55.11 | -55.10 | -55.07 | +| gjf full || -55.11 | -55.11 | -55.11 | -55.11 | -55.11 | -55.10 | -55.07 | +| langevin || -55.11 | -55.07 | -54.87 | -54.79 | -54.65 | X | X | +| nvt (Nose-Hoover) || -55.14 | -55.07 | -54.90 | -54.84 | -54.76 | X | X | +|-------------------||--------|--------|--------|--------|--------|--------|--------| + +KINETIC ENERGY (eV) +| Δt || 0.01 | 0.05 | 0.10 | 0.11 | 0.12 | 0.13 | 0.14 | +|===================||========|========|========|========|========|========|========| +| gjf half || 1.117 | 1.116 | 1.119 | 1.119 | 1.123 | 1.136 | 1.170 | +| gjf full || 1.116 | 1.071 | 0.938 | 0.898 | 0.858 | 0.817 | 0.780 | +| langevin || 1.110 | 1.113 | 1.121 | 1.129 | 1.157 | X | X | +| nvt (Nose-Hoover) || 1.083 | 1.109 | 1.112 | 1.113 | 1.114 | X | X | +|-------------------||--------|--------|--------|--------|--------|--------|--------| Script Commands: diff --git a/examples/gjf/in.gjf.vfull b/examples/gjf/in.gjf.vfull index 127b8fef87..c36a5b0690 100644 --- a/examples/gjf/in.gjf.vfull +++ b/examples/gjf/in.gjf.vfull @@ -17,9 +17,9 @@ compute myPE all pe fix lang all langevin/gjf 10 10 1 26488 vel vfull method 1 -thermo 2000 run 5000 -fix energies all ave/time 1 20000 20000 c_myKE c_myPE file ave.out +fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out +thermo 2000 run 35000 \ No newline at end of file diff --git a/examples/gjf/in.gjf.vhalf b/examples/gjf/in.gjf.vhalf index ba805d693a..0635998954 100644 --- a/examples/gjf/in.gjf.vhalf +++ b/examples/gjf/in.gjf.vhalf @@ -19,7 +19,7 @@ fix lang all langevin/gjf 10 10 1 26488 run 5000 -fix energies all ave/time 1 20000 20000 c_myKE c_myPE file ave.out +fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out thermo 2000 run 35000 \ No newline at end of file diff --git a/examples/gjf/log.2Apr25.gjf.vfull.g++.1 b/examples/gjf/log.2Apr25.gjf.vfull.g++.1 index 30a0f7e9c9..3bb204704b 100644 --- a/examples/gjf/log.2Apr25.gjf.vfull.g++.1 +++ b/examples/gjf/log.2Apr25.gjf.vfull.g++.1 @@ -1,4 +1,4 @@ -LAMMPS (2 Apr 2025 - Development - 09d3ac0a1b-modified) +LAMMPS (2 Apr 2025 - Development - e4c3b0c05e-modified) OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99) using 1 OpenMP thread(s) per MPI task # GJ thermostat @@ -55,7 +55,6 @@ compute myPE all pe fix lang all langevin/gjf 10 10 1 26488 vel vfull method 1 -thermo 2000 run 5000 CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE @@ -72,7 +71,7 @@ url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506}, doi = {10.1080/00268976.2019.1662506}, journal = {Molecular Physics}, author = {Grønbech-Jensen, Niels}, -year = {2020}, +year = {2020} } CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE @@ -93,24 +92,22 @@ Neighbor list info ... Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes Step Temp E_pair E_mol TotEng Press 0 10 -56.207652 0 -55.092137 33.341103 - 2000 8.5329244 -55.06102 0 -54.10916 343.01505 - 4000 7.8165794 -55.128803 0 -54.256851 322.07344 5000 8.4535562 -55.150518 0 -54.207511 318.20862 -Loop time of 2.26625 on 1 procs for 5000 steps with 864 atoms +Loop time of 2.32365 on 1 procs for 5000 steps with 864 atoms -Performance: 19062.331 ns/day, 0.001 hours/ns, 2206.288 timesteps/s, 1.906 Matom-step/s -100.0% CPU use with 1 MPI tasks x 1 OpenMP threads +Performance: 18591.401 ns/day, 0.001 hours/ns, 2151.782 timesteps/s, 1.859 Matom-step/s +99.9% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 1.2907 | 1.2907 | 1.2907 | 0.0 | 56.95 -Bond | 0.00032498 | 0.00032498 | 0.00032498 | 0.0 | 0.01 -Neigh | 0.28343 | 0.28343 | 0.28343 | 0.0 | 12.51 -Comm | 0.055491 | 0.055491 | 0.055491 | 0.0 | 2.45 -Output | 8.5614e-05 | 8.5614e-05 | 8.5614e-05 | 0.0 | 0.00 -Modify | 0.61129 | 0.61129 | 0.61129 | 0.0 | 26.97 -Other | | 0.02497 | | | 1.10 +Pair | 1.3067 | 1.3067 | 1.3067 | 0.0 | 56.23 +Bond | 0.00027472 | 0.00027472 | 0.00027472 | 0.0 | 0.01 +Neigh | 0.28218 | 0.28218 | 0.28218 | 0.0 | 12.14 +Comm | 0.057388 | 0.057388 | 0.057388 | 0.0 | 2.47 +Output | 5.6019e-05 | 5.6019e-05 | 5.6019e-05 | 0.0 | 0.00 +Modify | 0.65159 | 0.65159 | 0.65159 | 0.0 | 28.04 +Other | | 0.02549 | | | 1.10 Nlocal: 864 ave 864 max 864 min Histogram: 1 0 0 0 0 0 0 0 0 0 @@ -125,8 +122,9 @@ Ave special neighs/atom = 0 Neighbor list builds = 258 Dangerous builds = 0 -fix energies all ave/time 1 20000 20000 c_myKE c_myPE file ave.out +fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out +thermo 2000 run 35000 Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes @@ -150,21 +148,21 @@ Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes 36000 8.3048574 -55.079475 0 -54.153056 338.04291 38000 8.8708544 -55.108991 0 -54.119434 330.70097 40000 8.4012779 -55.080817 0 -54.143642 338.54326 -Loop time of 19.3814 on 1 procs for 35000 steps with 864 atoms +Loop time of 19.3834 on 1 procs for 35000 steps with 864 atoms -Performance: 15602.554 ns/day, 0.002 hours/ns, 1805.851 timesteps/s, 1.560 Matom-step/s -100.0% CPU use with 1 MPI tasks x 1 OpenMP threads +Performance: 15601.000 ns/day, 0.002 hours/ns, 1805.671 timesteps/s, 1.560 Matom-step/s +99.9% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 12.065 | 12.065 | 12.065 | 0.0 | 62.25 -Bond | 0.0021494 | 0.0021494 | 0.0021494 | 0.0 | 0.01 -Neigh | 2.3755 | 2.3755 | 2.3755 | 0.0 | 12.26 -Comm | 0.40916 | 0.40916 | 0.40916 | 0.0 | 2.11 -Output | 0.0004907 | 0.0004907 | 0.0004907 | 0.0 | 0.00 -Modify | 4.3528 | 4.3528 | 4.3528 | 0.0 | 22.46 -Other | | 0.1765 | | | 0.91 +Pair | 11.828 | 11.828 | 11.828 | 0.0 | 61.02 +Bond | 0.0024424 | 0.0024424 | 0.0024424 | 0.0 | 0.01 +Neigh | 2.3519 | 2.3519 | 2.3519 | 0.0 | 12.13 +Comm | 0.42251 | 0.42251 | 0.42251 | 0.0 | 2.18 +Output | 0.00056686 | 0.00056686 | 0.00056686 | 0.0 | 0.00 +Modify | 4.5987 | 4.5987 | 4.5987 | 0.0 | 23.73 +Other | | 0.1793 | | | 0.93 Nlocal: 864 ave 864 max 864 min Histogram: 1 0 0 0 0 0 0 0 0 0 diff --git a/examples/gjf/log.2Apr25.gjf.vfull.g++.4 b/examples/gjf/log.2Apr25.gjf.vfull.g++.4 index 454887a106..4c6fe82e54 100644 --- a/examples/gjf/log.2Apr25.gjf.vfull.g++.4 +++ b/examples/gjf/log.2Apr25.gjf.vfull.g++.4 @@ -1,4 +1,4 @@ -LAMMPS (2 Apr 2025 - Development - 09d3ac0a1b-modified) +LAMMPS (2 Apr 2025 - Development - e4c3b0c05e-modified) OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99) using 1 OpenMP thread(s) per MPI task # GJ thermostat @@ -55,7 +55,6 @@ compute myPE all pe fix lang all langevin/gjf 10 10 1 26488 vel vfull method 1 -thermo 2000 run 5000 CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE @@ -72,7 +71,7 @@ url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506}, doi = {10.1080/00268976.2019.1662506}, journal = {Molecular Physics}, author = {Grønbech-Jensen, Niels}, -year = {2020}, +year = {2020} } CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE @@ -93,24 +92,22 @@ Neighbor list info ... Per MPI rank memory allocation (min/avg/max) = 6.427 | 6.427 | 6.427 Mbytes Step Temp E_pair E_mol TotEng Press 0 10 -56.207652 0 -55.092137 33.341103 - 2000 8.63533 -55.154093 0 -54.190809 317.34773 - 4000 8.2482881 -55.091551 0 -54.171442 333.83766 5000 7.946377 -55.076514 0 -54.190084 337.31999 -Loop time of 2.06641 on 4 procs for 5000 steps with 864 atoms +Loop time of 2.0791 on 4 procs for 5000 steps with 864 atoms -Performance: 20905.811 ns/day, 0.001 hours/ns, 2419.654 timesteps/s, 2.091 Matom-step/s -62.1% CPU use with 4 MPI tasks x 1 OpenMP threads +Performance: 20778.210 ns/day, 0.001 hours/ns, 2404.885 timesteps/s, 2.078 Matom-step/s +64.9% CPU use with 4 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 0.53593 | 0.54393 | 0.55031 | 0.7 | 26.32 -Bond | 0.0004435 | 0.00052388 | 0.000586 | 0.0 | 0.03 -Neigh | 0.10844 | 0.11152 | 0.11687 | 1.0 | 5.40 -Comm | 0.94337 | 0.94742 | 0.94979 | 0.3 | 45.85 -Output | 0.00076423 | 0.00077499 | 0.00078949 | 0.0 | 0.04 -Modify | 0.24993 | 0.25116 | 0.25286 | 0.2 | 12.15 -Other | | 0.2111 | | | 10.21 +Pair | 0.52985 | 0.53809 | 0.54586 | 0.8 | 25.88 +Bond | 0.00054599 | 0.00057676 | 0.00061668 | 0.0 | 0.03 +Neigh | 0.10872 | 0.10965 | 0.11044 | 0.2 | 5.27 +Comm | 0.95581 | 0.96284 | 0.97025 | 0.5 | 46.31 +Output | 0.00027472 | 0.00028374 | 0.00029435 | 0.0 | 0.01 +Modify | 0.25077 | 0.25198 | 0.2529 | 0.2 | 12.12 +Other | | 0.2157 | | | 10.37 Nlocal: 216 ave 216 max 216 min Histogram: 4 0 0 0 0 0 0 0 0 0 @@ -125,8 +122,9 @@ Ave special neighs/atom = 0 Neighbor list builds = 273 Dangerous builds = 0 -fix energies all ave/time 1 20000 20000 c_myKE c_myPE file ave.out +fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out +thermo 2000 run 35000 Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule Per MPI rank memory allocation (min/avg/max) = 6.428 | 6.428 | 6.428 Mbytes @@ -150,21 +148,21 @@ Per MPI rank memory allocation (min/avg/max) = 6.428 | 6.428 | 6.428 Mbytes 36000 8.0936615 -55.152202 0 -54.249342 316.20169 38000 7.999652 -55.048407 0 -54.156034 345.07945 40000 8.6699753 -55.087634 0 -54.120485 337.23709 -Loop time of 17.4329 on 4 procs for 35000 steps with 864 atoms +Loop time of 17.5212 on 4 procs for 35000 steps with 864 atoms -Performance: 17346.528 ns/day, 0.001 hours/ns, 2007.700 timesteps/s, 1.735 Matom-step/s -65.5% CPU use with 4 MPI tasks x 1 OpenMP threads +Performance: 17259.111 ns/day, 0.001 hours/ns, 1997.582 timesteps/s, 1.726 Matom-step/s +63.6% CPU use with 4 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 5.1129 | 5.1423 | 5.1581 | 0.8 | 29.50 -Bond | 0.0029266 | 0.0034929 | 0.0038279 | 0.6 | 0.02 -Neigh | 0.85004 | 0.87543 | 0.90992 | 2.4 | 5.02 -Comm | 6.6875 | 6.6993 | 6.7123 | 0.4 | 38.43 -Output | 0.0043725 | 0.0045522 | 0.005017 | 0.4 | 0.03 -Modify | 3.2524 | 3.2579 | 3.2677 | 0.3 | 18.69 -Other | | 1.45 | | | 8.32 +Pair | 4.9924 | 5.0801 | 5.1446 | 2.4 | 28.99 +Bond | 0.0041975 | 0.0044234 | 0.0046863 | 0.3 | 0.03 +Neigh | 0.85805 | 0.86369 | 0.87021 | 0.5 | 4.93 +Comm | 6.7339 | 6.8008 | 6.8891 | 2.1 | 38.81 +Output | 0.0044786 | 0.0046306 | 0.0050572 | 0.4 | 0.03 +Modify | 3.2849 | 3.2919 | 3.2981 | 0.3 | 18.79 +Other | | 1.476 | | | 8.42 Nlocal: 216 ave 222 max 210 min Histogram: 2 0 0 0 0 0 0 0 0 2 diff --git a/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 b/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 index 5d07985f1c..7ea6355b64 100644 --- a/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 +++ b/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 @@ -1,4 +1,4 @@ -LAMMPS (2 Apr 2025 - Development - 09d3ac0a1b-modified) +LAMMPS (2 Apr 2025 - Development - e4c3b0c05e-modified) OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99) using 1 OpenMP thread(s) per MPI task # GJ thermostat @@ -21,7 +21,7 @@ Finding 1-2 1-3 1-4 neighbors ... 0 = max # of 1-4 neighbors 1 = max # of special neighbors special bonds CPU = 0.000 seconds - read_data CPU = 0.007 seconds + read_data CPU = 0.012 seconds include ff-argon.lmp ############################# #Atoms types - mass - charge# @@ -71,7 +71,20 @@ url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506}, doi = {10.1080/00268976.2019.1662506}, journal = {Molecular Physics}, author = {Grønbech-Jensen, Niels}, -year = {2020}, +year = {2020} +} + +- Langevin GJ-I vhalf method: doi:10.1080/00268976.2019.1570369 + +@Article{jensen_accurate_2019, +title = {Accurate configurational and kinetic statistics in discrete-time Langevin systems}, +volume = {117}, +url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1570369}, +doi = {10.1080/00268976.2019.1570369}, +number = {18}, +journal = {Molecular Physics}, +author = {Jensen, Lucas Frese Grønbech and Grønbech-Jensen, Niels}, +year = {2019} } CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE @@ -93,21 +106,21 @@ Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes Step Temp E_pair E_mol TotEng Press 0 10 -56.207652 0 -55.092137 33.341103 5000 9.7731898 -55.150518 0 -54.060304 322.94195 -Loop time of 2.32188 on 1 procs for 5000 steps with 864 atoms +Loop time of 2.33802 on 1 procs for 5000 steps with 864 atoms -Performance: 18605.621 ns/day, 0.001 hours/ns, 2153.428 timesteps/s, 1.861 Matom-step/s -99.8% CPU use with 1 MPI tasks x 1 OpenMP threads +Performance: 18477.199 ns/day, 0.001 hours/ns, 2138.565 timesteps/s, 1.848 Matom-step/s +99.6% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 1.3154 | 1.3154 | 1.3154 | 0.0 | 56.65 -Bond | 0.00030719 | 0.00030719 | 0.00030719 | 0.0 | 0.01 -Neigh | 0.28939 | 0.28939 | 0.28939 | 0.0 | 12.46 -Comm | 0.055279 | 0.055279 | 0.055279 | 0.0 | 2.38 -Output | 4.3714e-05 | 4.3714e-05 | 4.3714e-05 | 0.0 | 0.00 -Modify | 0.63675 | 0.63675 | 0.63675 | 0.0 | 27.42 -Other | | 0.02471 | | | 1.06 +Pair | 1.3074 | 1.3074 | 1.3074 | 0.0 | 55.92 +Bond | 0.00028207 | 0.00028207 | 0.00028207 | 0.0 | 0.01 +Neigh | 0.27625 | 0.27625 | 0.27625 | 0.0 | 11.82 +Comm | 0.056554 | 0.056554 | 0.056554 | 0.0 | 2.42 +Output | 5.3527e-05 | 5.3527e-05 | 5.3527e-05 | 0.0 | 0.00 +Modify | 0.67274 | 0.67274 | 0.67274 | 0.0 | 28.77 +Other | | 0.02477 | | | 1.06 Nlocal: 864 ave 864 max 864 min Histogram: 1 0 0 0 0 0 0 0 0 0 @@ -122,7 +135,7 @@ Ave special neighs/atom = 0 Neighbor list builds = 258 Dangerous builds = 0 -fix energies all ave/time 1 20000 20000 c_myKE c_myPE file ave.out +fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out thermo 2000 run 35000 @@ -148,21 +161,21 @@ Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes 36000 9.7252627 -55.079475 0 -53.994608 343.13769 38000 10.438984 -55.108991 0 -53.944506 336.32562 40000 10.238268 -55.080817 0 -53.938723 345.13228 -Loop time of 19.5673 on 1 procs for 35000 steps with 864 atoms +Loop time of 19.4518 on 1 procs for 35000 steps with 864 atoms -Performance: 15454.339 ns/day, 0.002 hours/ns, 1788.697 timesteps/s, 1.545 Matom-step/s +Performance: 15546.156 ns/day, 0.002 hours/ns, 1799.324 timesteps/s, 1.555 Matom-step/s 99.9% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 12.065 | 12.065 | 12.065 | 0.0 | 61.66 -Bond | 0.002157 | 0.002157 | 0.002157 | 0.0 | 0.01 -Neigh | 2.384 | 2.384 | 2.384 | 0.0 | 12.18 -Comm | 0.40912 | 0.40912 | 0.40912 | 0.0 | 2.09 -Output | 0.00050573 | 0.00050573 | 0.00050573 | 0.0 | 0.00 -Modify | 4.5321 | 4.5321 | 4.5321 | 0.0 | 23.16 -Other | | 0.1749 | | | 0.89 +Pair | 11.789 | 11.789 | 11.789 | 0.0 | 60.61 +Bond | 0.0022933 | 0.0022933 | 0.0022933 | 0.0 | 0.01 +Neigh | 2.3359 | 2.3359 | 2.3359 | 0.0 | 12.01 +Comm | 0.4146 | 0.4146 | 0.4146 | 0.0 | 2.13 +Output | 0.0005505 | 0.0005505 | 0.0005505 | 0.0 | 0.00 +Modify | 4.7341 | 4.7341 | 4.7341 | 0.0 | 24.34 +Other | | 0.1749 | | | 0.90 Nlocal: 864 ave 864 max 864 min Histogram: 1 0 0 0 0 0 0 0 0 0 diff --git a/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 b/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 index cac135e578..088fa73044 100644 --- a/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 +++ b/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 @@ -1,4 +1,4 @@ -LAMMPS (2 Apr 2025 - Development - 09d3ac0a1b-modified) +LAMMPS (2 Apr 2025 - Development - e4c3b0c05e-modified) OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99) using 1 OpenMP thread(s) per MPI task # GJ thermostat @@ -21,7 +21,7 @@ Finding 1-2 1-3 1-4 neighbors ... 0 = max # of 1-4 neighbors 1 = max # of special neighbors special bonds CPU = 0.002 seconds - read_data CPU = 0.014 seconds + read_data CPU = 0.015 seconds include ff-argon.lmp ############################# #Atoms types - mass - charge# @@ -71,7 +71,20 @@ url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506}, doi = {10.1080/00268976.2019.1662506}, journal = {Molecular Physics}, author = {Grønbech-Jensen, Niels}, -year = {2020}, +year = {2020} +} + +- Langevin GJ-I vhalf method: doi:10.1080/00268976.2019.1570369 + +@Article{jensen_accurate_2019, +title = {Accurate configurational and kinetic statistics in discrete-time Langevin systems}, +volume = {117}, +url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1570369}, +doi = {10.1080/00268976.2019.1570369}, +number = {18}, +journal = {Molecular Physics}, +author = {Jensen, Lucas Frese Grønbech and Grønbech-Jensen, Niels}, +year = {2019} } CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE @@ -93,21 +106,21 @@ Per MPI rank memory allocation (min/avg/max) = 6.427 | 6.427 | 6.427 Mbytes Step Temp E_pair E_mol TotEng Press 0 10 -56.207652 0 -55.092137 33.341103 5000 9.3726166 -55.076514 0 -54.030985 342.43571 -Loop time of 2.09273 on 4 procs for 5000 steps with 864 atoms +Loop time of 2.09482 on 4 procs for 5000 steps with 864 atoms -Performance: 20642.908 ns/day, 0.001 hours/ns, 2389.226 timesteps/s, 2.064 Matom-step/s -62.7% CPU use with 4 MPI tasks x 1 OpenMP threads +Performance: 20622.269 ns/day, 0.001 hours/ns, 2386.837 timesteps/s, 2.062 Matom-step/s +64.9% CPU use with 4 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 0.53092 | 0.54173 | 0.55017 | 1.0 | 25.89 -Bond | 0.00049283 | 0.00052862 | 0.00057217 | 0.0 | 0.03 -Neigh | 0.10693 | 0.11084 | 0.1158 | 1.0 | 5.30 -Comm | 0.95622 | 0.96088 | 0.96553 | 0.4 | 45.92 -Output | 0.00027302 | 0.00027645 | 0.00028114 | 0.0 | 0.01 -Modify | 0.2615 | 0.26289 | 0.2635 | 0.2 | 12.56 -Other | | 0.2156 | | | 10.30 +Pair | 0.5344 | 0.54263 | 0.54877 | 0.8 | 25.90 +Bond | 0.00049036 | 0.00055513 | 0.00061267 | 0.0 | 0.03 +Neigh | 0.10741 | 0.10898 | 0.11004 | 0.3 | 5.20 +Comm | 0.95869 | 0.96388 | 0.96936 | 0.4 | 46.01 +Output | 0.00025982 | 0.00026553 | 0.00026761 | 0.0 | 0.01 +Modify | 0.26091 | 0.26108 | 0.26137 | 0.0 | 12.46 +Other | | 0.2174 | | | 10.38 Nlocal: 216 ave 216 max 216 min Histogram: 4 0 0 0 0 0 0 0 0 0 @@ -122,7 +135,7 @@ Ave special neighs/atom = 0 Neighbor list builds = 273 Dangerous builds = 0 -fix energies all ave/time 1 20000 20000 c_myKE c_myPE file ave.out +fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out thermo 2000 run 35000 @@ -148,21 +161,21 @@ Per MPI rank memory allocation (min/avg/max) = 6.428 | 6.428 | 6.428 Mbytes 36000 9.6032778 -55.152202 0 -54.080942 321.61646 38000 9.8802995 -55.048407 0 -53.946245 351.82506 40000 10.372288 -55.087634 0 -53.93059 343.34304 -Loop time of 17.6209 on 4 procs for 35000 steps with 864 atoms +Loop time of 17.6294 on 4 procs for 35000 steps with 864 atoms -Performance: 17161.415 ns/day, 0.001 hours/ns, 1986.275 timesteps/s, 1.716 Matom-step/s -64.2% CPU use with 4 MPI tasks x 1 OpenMP threads +Performance: 17153.172 ns/day, 0.001 hours/ns, 1985.321 timesteps/s, 1.715 Matom-step/s +63.7% CPU use with 4 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 5.0439 | 5.1029 | 5.1534 | 2.0 | 28.96 -Bond | 0.0031373 | 0.0034925 | 0.0038794 | 0.5 | 0.02 -Neigh | 0.84082 | 0.87558 | 0.91665 | 3.0 | 4.97 -Comm | 6.7337 | 6.7936 | 6.8487 | 1.6 | 38.55 -Output | 0.0044073 | 0.004559 | 0.0049773 | 0.4 | 0.03 -Modify | 3.3556 | 3.3647 | 3.3722 | 0.3 | 19.10 -Other | | 1.476 | | | 8.38 +Pair | 5.0222 | 5.1025 | 5.1759 | 2.5 | 28.94 +Bond | 0.0040066 | 0.0042899 | 0.0045178 | 0.3 | 0.02 +Neigh | 0.8511 | 0.85865 | 0.86632 | 0.6 | 4.87 +Comm | 6.7448 | 6.8089 | 6.8793 | 1.8 | 38.62 +Output | 0.0043966 | 0.0045518 | 0.0049846 | 0.4 | 0.03 +Modify | 3.364 | 3.3672 | 3.3717 | 0.2 | 19.10 +Other | | 1.483 | | | 8.41 Nlocal: 216 ave 222 max 210 min Histogram: 2 0 0 0 0 0 0 0 0 2 diff --git a/src/EXTRA-FIX/fix_langevin_gjf.cpp b/src/EXTRA-FIX/fix_langevin_gjf.cpp index ea0e75c598..73ab2a852f 100644 --- a/src/EXTRA-FIX/fix_langevin_gjf.cpp +++ b/src/EXTRA-FIX/fix_langevin_gjf.cpp @@ -41,7 +41,6 @@ using namespace FixConst; enum { NOBIAS, BIAS }; enum { CONSTANT, EQUAL, ATOM }; - static const char cite_langevin_gjf[] = "Langevin GJ methods: doi:10.1080/00268976.2019.1662506\n\n" "@Article{gronbech-jensen_complete_2020,\n" @@ -52,7 +51,7 @@ static const char cite_langevin_gjf[] = "doi = {10.1080/00268976.2019.1662506},\n" "journal = {Molecular Physics},\n" "author = {Grønbech-Jensen, Niels},\n" - "year = {2020},\n" + "year = {2020}\n" "}\n\n"; static const char cite_langevin_gjf_7[] = @@ -63,7 +62,6 @@ static const char cite_langevin_gjf_7[] = "number = {18},\n" "url = {https://doi.org/10.1063/5.0066008},\n" "doi = {10.1063/5.0066008},\n" - "urldate = {2021-11-14},\n" "journal = {J. Chem. Phys.},\n" "author = {Finkelstein, Joshua and Cheng, Chungho and Fiorin, Giacomo and Seibold, Benjamin and Grønbech-Jensen, Niels},\n" "year = {2021},\n" @@ -78,12 +76,24 @@ static const char cite_langevin_gjf_8[] = "number = {10},\n" "url = {https://doi.org/10.1007/s10955-024-03345-1},\n" "doi = {10.1007/s10955-024-03345-1},\n" - "urldate = {2024-10-22},\n" "journal = {J. Stat. Phys.},\n" "author = {Gronbech-Jensen, Niels},\n" "year = {2024},\n" "pages = {137}\n" "}\n\n"; + +static const char cite_langevin_gjf_vhalf[] = + "Langevin GJ-I vhalf method: doi:10.1080/00268976.2019.1570369\n\n" + "@Article{jensen_accurate_2019,\n" + "title = {Accurate configurational and kinetic statistics in discrete-time Langevin systems},\n" + "volume = {117},\n" + "url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1570369},\n" + "doi = {10.1080/00268976.2019.1570369},\n" + "number = {18},\n" + "journal = {Molecular Physics},\n" + "author = {Jensen, Lucas Frese Grønbech and Grønbech-Jensen, Niels},\n" + "year = {2019}\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ @@ -155,6 +165,7 @@ FixLangevinGJF::FixLangevinGJF(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR, "Illegal fix langevin/gjf command"); } + if (GJmethod == 1 && osflag == 0) if (lmp->citeme) lmp->citeme->add(cite_langevin_gjf_vhalf); // set temperature = nullptr, user can override via fix_modify if wants bias id_temp = nullptr; @@ -236,6 +247,11 @@ void FixLangevinGJF::init() error->all(FLERR, "Fix langevin/gjf and run style respa are not compatible"); } + if (temperature && temperature->tempbias) + tbiasflag = BIAS; + else + tbiasflag = NOBIAS; + // Complete set of thermostats is given in Gronbech-Jensen, Molecular Physics, 118 (2020) switch (GJmethod) { case 1: @@ -320,58 +336,95 @@ void FixLangevinGJF::initial_integrate(int /* vflag */) } compute_target(); + if (tbiasflag) temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - - if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); - if (rmass) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); m = rmass[i]; beta = tsqrt * sqrt(2.0*dt*m*boltz/t_period/mvv2e) / ftm2v; - } else { + + fran[0] = beta*random->gaussian(); + fran[1] = beta*random->gaussian(); + fran[2] = beta*random->gaussian(); + + // First integration delivers Eq. 24a and 24b: + dtfm = dtf / m; + v[i][0] += csq * dtfm * f[i][0]; + v[i][1] += csq * dtfm * f[i][1]; + v[i][2] += csq * dtfm * f[i][2]; + x[i][0] += 0.5 * csq * dt * v[i][0]; + x[i][1] += 0.5 * csq * dt * v[i][1]; + x[i][2] += 0.5 * csq * dt * v[i][2]; + + if (tbiasflag) temperature->remove_bias(i, v[i]); + + // Calculate Eq. 24c: + lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c3sqrt / (2.0 * m)) * fran[0]; + lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c3sqrt / (2.0 * m)) * fran[1]; + lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c3sqrt / (2.0 * m)) * fran[2]; + + // Calculate Eq. 24d + v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * csq * (0.5 / m) * fran[0]; + v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * csq * (0.5 / m) * fran[1]; + v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * csq * (0.5 / m) * fran[2]; + + if (tbiasflag) temperature->restore_bias(i, v[i]); + if (tbiasflag) temperature->restore_bias(i, lv[i]); + + // Calculate Eq. 24e. Final integrator then calculates Eq. 24f after force update. + x[i][0] += 0.5 * csq * dt * v[i][0]; + x[i][1] += 0.5 * csq * dt * v[i][1]; + x[i][2] += 0.5 * csq * dt * v[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); m = mass[type[i]]; beta = tsqrt * sqrt(2.0*dt*m*boltz/t_period/mvv2e) / ftm2v; + + fran[0] = beta*random->gaussian(); + fran[1] = beta*random->gaussian(); + fran[2] = beta*random->gaussian(); + + // First integration delivers Eq. 24a and 24b: + dtfm = dtf / m; + v[i][0] += csq * dtfm * f[i][0]; + v[i][1] += csq * dtfm * f[i][1]; + v[i][2] += csq * dtfm * f[i][2]; + x[i][0] += 0.5 * csq * dt * v[i][0]; + x[i][1] += 0.5 * csq * dt * v[i][1]; + x[i][2] += 0.5 * csq * dt * v[i][2]; + + if (tbiasflag) temperature->remove_bias(i, v[i]); + + // Calculate Eq. 24c: + lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c3sqrt / (2.0 * m)) * fran[0]; + lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c3sqrt / (2.0 * m)) * fran[1]; + lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c3sqrt / (2.0 * m)) * fran[2]; + + // Calculate Eq. 24d + v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * csq * (0.5 / m) * fran[0]; + v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * csq * (0.5 / m) * fran[1]; + v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * csq * (0.5 / m) * fran[2]; + + if (tbiasflag) temperature->restore_bias(i, v[i]); + if (tbiasflag) temperature->restore_bias(i, lv[i]); + + // Calculate Eq. 24e. Final integrator then calculates Eq. 24f after force update. + x[i][0] += 0.5 * csq * dt * v[i][0]; + x[i][1] += 0.5 * csq * dt * v[i][1]; + x[i][2] += 0.5 * csq * dt * v[i][2]; } - - fran[0] = beta*random->gaussian(); - fran[1] = beta*random->gaussian(); - fran[2] = beta*random->gaussian(); - - // First integration delivers Eq. 24a and 24b: - dtfm = dtf / m; - v[i][0] += csq * dtfm * f[i][0]; - v[i][1] += csq * dtfm * f[i][1]; - v[i][2] += csq * dtfm * f[i][2]; - x[i][0] += 0.5 * csq * dt * v[i][0]; - x[i][1] += 0.5 * csq * dt * v[i][1]; - x[i][2] += 0.5 * csq * dt * v[i][2]; - - // Calculate Eq. 24c: - lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c3sqrt / (2.0 * m)) * fran[0]; - lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c3sqrt / (2.0 * m)) * fran[1]; - lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c3sqrt / (2.0 * m)) * fran[2]; - - // Calculate Eq. 24d - v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * csq * (0.5 / m) * fran[0]; - v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * csq * (0.5 / m) * fran[1]; - v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * csq * (0.5 / m) * fran[2]; - - // Calculate Eq. 24e. Final integrator then calculates Eq. 24f after force update. - x[i][0] += 0.5 * csq * dt * v[i][0]; - x[i][1] += 0.5 * csq * dt * v[i][1]; - x[i][2] += 0.5 * csq * dt * v[i][2]; } } } void FixLangevinGJF::final_integrate() { - double dtfm; - double dt = update->dt; - double ftm2v = force->ftm2v; - double dtf = 0.5 * dt * ftm2v; - double csq = sqrt(gjfc3 / gjfc1); - double **v = atom->v; double **f = atom->f; double *rmass = atom->rmass; @@ -381,6 +434,10 @@ void FixLangevinGJF::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; + double dtfm; + double dtf = 0.5 * update->dt * force->ftm2v; + double csq = sqrt(gjfc3 / gjfc1); + // Calculate Eq. 24f. if (rmass) { for (int i = 0; i < nlocal; i++) diff --git a/src/EXTRA-FIX/fix_langevin_gjf.h b/src/EXTRA-FIX/fix_langevin_gjf.h index 4b65051cf5..45ccf4cbe2 100644 --- a/src/EXTRA-FIX/fix_langevin_gjf.h +++ b/src/EXTRA-FIX/fix_langevin_gjf.h @@ -44,7 +44,7 @@ class FixLangevinGJF : public Fix { int unpack_exchange(int, double *) override; protected: - int osflag, GJmethod, maxatom, lv_allocated; + int osflag, tbiasflag, GJmethod, maxatom, lv_allocated; double t_start, t_stop, t_period, t_target, tsqrt; double gjfc1, gjfc2, gjfc3; int tstyle, tvar;