Rethinking thermostats

examples/Argon/Argon.cpp

@@ -1,35 +1,33 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-#include <format>
-
 #include <CLI/App.hpp>
 #include <CLI/Config.hpp>
 #include <CLI/Formatter.hpp>
 #include <Kokkos_Core.hpp>
 #include <algorithm>
+#include <format>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 
 #include "Cabana_NeighborList.hpp"
-#include "action/BerendsenThermostat.hpp"
-#include "action/LangevinThermostat.hpp"
 #include "action/LennardJones.hpp"
 #include "action/LimitAcceleration.hpp"
 #include "action/LimitVelocity.hpp"
 #include "action/VelocityVerlet.hpp"
+#include "action/VelocityVerletLangevinThermostat.hpp"
 #include "analysis/KineticEnergy.hpp"
 #include "analysis/MeanSquareDisplacement.hpp"
 #include "analysis/Pressure.hpp"
@@ -72,7 +70,7 @@ struct Config
 
     idx_t estimatedMaxNeighbors = 60;
 
-    const std::string resName = "Argon"; 
+    const std::string resName = "Argon";
     const std::vector<std::string> typeNames = {"Ar"};
 };
 
@@ -123,11 +121,18 @@ void LJ(Config& config)
     auto rho = real_c(atoms.numLocalAtoms) / volume;
     std::cout << "rho: " << rho << std::endl;
 
-    io::dumpGRO("atoms_initial.gro", atoms, subdomain, 0_r, "Argon", config.resName, config.typeNames, false);
+    io::dumpGRO("atoms_initial.gro",
+                atoms,
+                subdomain,
+                0_r,
+                "Argon",
+                config.resName,
+                config.typeNames,
+                false);
 
     communication::GhostLayer ghostLayer;
     action::LennardJones LJ(config.rc, config.sigma, config.epsilon, 0.7_r * config.sigma);
-    action::LangevinThermostat langevinThermostat(config.gamma, config.temperature, config.dt);
+    action::VelocityVerletLangevinThermostat langevinIntegrator(config.gamma, config.temperature);
     analysis::MeanSquareDisplacement meanSquareDisplacement;
     meanSquareDisplacement.reset(atoms);
     auto msd = 0_r;
@@ -141,7 +146,14 @@ void LJ(Config& config)
     std::ofstream fStat("statistics.txt");
     for (auto step = 0; step < config.nsteps; ++step)
     {
-        maxAtomDisplacement += action::VelocityVerlet::preForceIntegrate(atoms, config.dt);
+        if (step % 1000 != 0)
+        {
+            maxAtomDisplacement += action::VelocityVerlet::preForceIntegrate(atoms, config.dt);
+        }
+        else
+        {
+            maxAtomDisplacement += langevinIntegrator.preForceIntegrate(atoms, config.dt);
+        }
 
         if (maxAtomDisplacement >= config.skin * 0.5_r)
         {
@@ -222,21 +234,6 @@ void LJ(Config& config)
             //            false);
         }
 
-        if (step % 1000 == 0)
-        {
-            msd = meanSquareDisplacement.calc(atoms, subdomain) / (1000_r * config.dt);
-            if ((config.temperature > 0_r) && (step > 5000))
-            {
-                config.temperature -= 7.8e-3_r;
-                config.temperature = std::max(config.temperature, 0_r);
-            }
-
-            langevinThermostat.set(config.gamma, config.temperature * 0.5_r, config.dt);
-            langevinThermostat.apply(atoms);
-
-            meanSquareDisplacement.reset(atoms);
-        }
-
         action::VelocityVerlet::postForceIntegrate(atoms, config.dt);
     }
     fStat.close();

examples/Argon/MultiResArgon.cpp

@@ -1,37 +1,36 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-#include <format>
-
 #include <CLI/App.hpp>
 #include <CLI/Config.hpp>
 #include <CLI/Formatter.hpp>
 #include <Kokkos_Core.hpp>
 #include <algorithm>
+#include <format>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 
 #include "Cabana_NeighborList.hpp"
 #include "action/BerendsenThermostat.hpp"
 #include "action/ContributeMoleculeForceToAtoms.hpp"
-#include "action/LangevinThermostat.hpp"
 #include "action/LennardJones.hpp"
 #include "action/LimitAcceleration.hpp"
 #include "action/LimitVelocity.hpp"
 #include "action/UpdateMolecules.hpp"
 #include "action/VelocityVerlet.hpp"
+#include "action/VelocityVerletLangevinThermostat.hpp"
 #include "analysis/KineticEnergy.hpp"
 #include "analysis/MeanSquareDisplacement.hpp"
 #include "analysis/Pressure.hpp"
@@ -76,9 +75,8 @@ struct Config
 
     idx_t estimatedMaxNeighbors = 60;
 
-    const std::string resName = "Argon"; 
+    const std::string resName = "Argon";
     const std::vector<std::string> typeNames = {"Ar"};
-
 };
 
 data::Atoms fillDomainWithAtomsSC(const data::Subdomain& subdomain,
@@ -131,12 +129,19 @@ void LJ(Config& config)
     auto rho = real_c(atoms.numLocalAtoms) / volume;
     std::cout << "rho: " << rho << std::endl;
 
-    io::dumpGRO("atoms_initial.gro", atoms, subdomain, 0_r, "Argon", config.resName, config.typeNames, false);
+    io::dumpGRO("atoms_initial.gro",
+                atoms,
+                subdomain,
+                0_r,
+                "Argon",
+                config.resName,
+                config.typeNames,
+                false);
 
     communication::MultiResGhostLayer ghostLayer;
     weighting_function::Slab weightingFunction({-100_r, -100_r, -100_r}, 1_r, 1_r, 1);
     action::LennardJones LJ(config.rc, config.sigma, config.epsilon, 0.7_r * config.sigma);
-    action::LangevinThermostat langevinThermostat(config.gamma, config.temperature, config.dt);
+    action::VelocityVerletLangevinThermostat langevinIntegrator(config.gamma, config.temperature);
     analysis::MeanSquareDisplacement meanSquareDisplacement;
     meanSquareDisplacement.reset(atoms);
     auto msd = 0_r;
@@ -150,7 +155,14 @@ void LJ(Config& config)
     std::ofstream fStat("statistics.txt");
     for (auto step = 0; step < config.nsteps; ++step)
     {
-        maxAtomDisplacement += action::VelocityVerlet::preForceIntegrate(atoms, config.dt);
+        if (step % 1000 != 0)
+        {
+            maxAtomDisplacement += action::VelocityVerlet::preForceIntegrate(atoms, config.dt);
+        }
+        else
+        {
+            maxAtomDisplacement += langevinIntegrator.preForceIntegrate(atoms, config.dt);
+        }
 
         action::UpdateMolecules::update(molecules, atoms, weightingFunction);
 
@@ -230,29 +242,14 @@ void LJ(Config& config)
                         subdomain,
                         step * config.dt,
                         "Argon",
-                        config.resName, 
+                        config.resName,
                         config.typeNames,
                         false);
 
             //            io::dumpCSV("atoms_" + std::to_string(step) + ".csv", atoms,
             //            false);
         }
 
-        if (step % 1000 == 0)
-        {
-            msd = meanSquareDisplacement.calc(atoms, subdomain) / (1000_r * config.dt);
-            if ((config.temperature > 0_r) && (step > 5000))
-            {
-                config.temperature -= 7.8e-3_r;
-                config.temperature = std::max(config.temperature, 0_r);
-            }
-
-            langevinThermostat.set(config.gamma, config.temperature * 0.5_r, config.dt);
-            langevinThermostat.apply(atoms);
-
-            meanSquareDisplacement.reset(atoms);
-        }
-
         action::VelocityVerlet::postForceIntegrate(atoms, config.dt);
     }
     fStat.close();

examples/BinaryLennardJones/BinaryLennardJones.cpp

@@ -1,11 +1,11 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -72,7 +72,7 @@ int main(int argc, char* argv[])  // NOLINT
     std::cout << "==================" << std::endl;
 
     LJ(cfg);
-    
+
     MPI_Finalize();
 
     return EXIT_SUCCESS;

examples/BinaryLennardJones/NPT.cpp

@@ -1,11 +1,11 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -14,9 +14,8 @@
 
 #include "NPT.hpp"
 
-#include <format>
-
 #include <algorithm>
+#include <format>
 
 #include "action/BerendsenBarostat.hpp"
 #include "action/BerendsenThermostat.hpp"

examples/BinaryLennardJones/NPT.hpp

@@ -1,11 +1,11 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.

examples/BinaryLennardJones/NVT.cpp

@@ -1,11 +1,11 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -14,9 +14,8 @@
 
 #include "NVT.hpp"
 
-#include <format>
-
 #include <algorithm>
+#include <format>
 
 #include "action/BerendsenThermostat.hpp"
 #include "action/LennardJones.hpp"

examples/BinaryLennardJones/NVT.hpp

@@ -1,11 +1,11 @@
 // Copyright 2024 Sebastian Eibl
-// 
+//
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
-// 
+//
 //     https://www.apache.org/licenses/LICENSE-2.0
-// 
+//
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.