Chunglu development merge

.gitignore

@@ -12,4 +12,9 @@ __pycache__
 Testing/Temporary/CTestCostData.txt
 .eggs
 wheelhouse
-vcpkg_installed
+vcpkg_installed
+
+result
+build_*
+result
+result-*

.vscode/launch.json

@@ -4,6 +4,14 @@
   // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
   "version": "0.2.0",
   "configurations": [
+    {
+      "name": "Python Debugger: Current File",
+      "type": "debugpy",
+      "request": "launch",
+      "program": "${file}",
+      "console": "integratedTerminal",
+      "cwd": "${workspaceFolder}/scripts"
+    },
     {
       "name": "(ctest) Launch",
       "type": "cppdbg",

CMakeLists.txt

@@ -101,14 +101,18 @@ endif()
 ######### Dependencies
 ######################################################################
 
-set(Boost_USE_STATIC_LIBS ON)
 set(Boost_USE_MULTITHREADED ON)
 set(Boost_USE_STATIC_RUNTIME OFF)
-find_package(Boost 1.66.0 COMPONENTS system filesystem program_options unit_test_framework)
+
+# Try static Boost first, then fallback to shared Boost
+find_package(Boost 1.66.0 COMPONENTS filesystem program_options unit_test_framework)
 if(NOT Boost_FOUND)
-  message(WARNING "Static Boost not found, trying shared Boost")
   set(Boost_USE_STATIC_LIBS OFF)
-  find_package(Boost 1.66.0 REQUIRED COMPONENTS system filesystem program_options unit_test_framework)
+  find_package(Boost 1.66.0 REQUIRED COMPONENTS filesystem program_options unit_test_framework)
+endif()
+
+if(NOT TARGET Boost::system)
+  add_library(Boost::system INTERFACE IMPORTED)
 endif()
 
 add_subdirectory(src/magnet)

derivation.nix

@@ -2,7 +2,9 @@
 # To install `nix-env -u -f default.nix`
 # To develop `nix-shell` (will build the shell with dependencies)
 # To test build `nix-build`
-{ pkgs, python3 }:
+{ pkgs, python3,
+  visualiser ? false,
+  }:
 python3.pkgs.buildPythonPackage rec {
   name = "pydynamo";
   src = ./.;
@@ -29,15 +31,17 @@ python3.pkgs.buildPythonPackage rec {
     gcc
     pkg-config
     clang-tools
+  ]  ++ propagatedBuildInputs
+  ++ (lib.optionals visualiser [
     wrapGAppsHook3
-  ]  ++ propagatedBuildInputs;
+  ]);
 
   buildInputs = with pkgs; [
     # Basic build dependencies
-    bzip2.dev
-    boost.dev
+    bzip2
+    boost
     eigen
-    
+  ] ++ (lib.optionals visualiser [
     # Visualiser
     libGL
     gtkmm3.dev
@@ -47,5 +51,5 @@ python3.pkgs.buildPythonPackage rec {
     cairomm.dev
     libpng
     mesa
-  ];
+  ]);
 }

flake.nix

@@ -0,0 +1,26 @@
+{
+  description = "PyDynamO / DynamO build environment";
+
+  inputs = {
+    nixpkgs.url = "github:NixOS/nixpkgs/nixos-24.11";
+    flake-utils.url = "github:numtide/flake-utils";
+  };
+
+  outputs = { self, nixpkgs, flake-utils }:
+    flake-utils.lib.eachDefaultSystem (system:
+      let
+        pkgs = import nixpkgs {
+          inherit system;
+        };
+        pydynamo = pkgs.callPackage ./derivation.nix {};
+      in
+      {
+        packages.default = pydynamo;
+        packages.pydynamo = pydynamo;
+
+        devShells.default = pkgs.mkShell {
+          inputsFrom = [ pydynamo ];
+        };
+      }
+    );
+}

scripts/HS_stats.py

@@ -0,0 +1,155 @@
+#!/usr/bin/env python3
+import pydynamo
+from pydynamo import ET
+
+
+def setup_worker( config,                # The name of the config file to generate.
+                  state,                 # A dictionary of state variables to use
+                  logfile,               # File handle where to write progress/logging output
+                  particle_equil_events, # How many events will be run per particle to equilibrate the config. Useful if in setup you also need to equilibrate an intermediate configuration.
+):
+    from subprocess import check_call
+
+    #Here we work out how many unit cells to make the system out of for various packings
+    state = dict(state)
+    if 'InitState' not in state:
+        state['InitState'] = "FCC"
+
+    unitcellN = {
+        "FCC":4,
+        "HCP":4,
+        "BCC":2,
+        "SC":1,
+    }
+    Ncells_unrounded = (state['N'] / unitcellN[state['InitState']]) ** (1.0 / 3.0)
+    Ncells = int(round(Ncells_unrounded))
+    if abs(Ncells - Ncells_unrounded) > 0.1:
+        raise RuntimeError("Could not make "+str(state['N'])+" particles in an "+state['InitState']+" packing")
+
+    # Here, for tethered systems, we do not simulate state points if
+    # its going to be boring and "ideal". I only have worked out the
+    # spacing expression for FCC, so all other crystals will just be
+    # run regardless
+
+    if ("Rso" in state) and (state['Rso'] != float('inf')) and (state['InitState'] == "FCC"):
+        effrho = state['ndensity']*(state['Lambda']**3)
+        minR = max(0, (2**(2.5)*effrho)**(-1/3.0) - 0.5)
+        #phiT= state['ndensity'] * (4/3) * math.pi * minR**3
+        #minRho = max(0, (2**(1/6.0)-(6*state['ndensity']*(4/3)*minR**3)**(1/3))**3)
+        if state['Rso'] <= minR:
+            raise pydynamo.SkipThisPoint()
+
+    # This check is halting systems deep in the solid region, which should not be done!
+    #
+    #
+    ### Again, for tethered systems in FCC lattices we do not simulate
+    ### much beyond a multiple of the minimum tether radius.
+    ##if ("Rso" in state) and (state['ndensity'] >= 1.0) and (state['InitState'] == "FCC"):
+    ##    minR = max(0, (2**(2.5)*state['ndensity'])**(-1/3.0) - 0.5)
+    ##    if state['Rso'] >= 10*minR:
+    ##        raise pydynamo.SkipThisPoint()
+
+
+    # Thermostat
+    options = ''
+    if 'kT' in state:
+        if state['kT'] != float('inf'):
+            options = options + ' -T '+repr(state['kT'])
+        else:#infite temperature is a special case, we set well energies to zero
+            options = options + ' -T 1.0'
+
+    # Crystal lattice packing
+    packmode = {
+        'FCC':0,
+        'BCC':1,
+        'SC': 2,
+        'HCP':3,
+    }
+    options = options + ' --i1 '+str(packmode[state['InitState']]) +' -C '+str(Ncells)
+
+    # Square well or hard sphere?
+    if state['Lambda'] != float('inf'):
+        if state['kT'] != float('inf'):
+            options = options + ' -m 1 --f1 '+repr(state['Lambda'])
+        else: # infinite temperature is a special case, we set well energies to zero
+            options = options + ' -m 1 --f1 '+repr(state['Lambda']) + " --f2 0.0"
+    else:
+        options = options + ' -m 0'
+
+    # denisty
+    options = options + ' -d ' + repr(state['ndensity'])
+
+    # Execution of dynamod
+    print('# dynamod'+options+' -o '+config, file=logfile)
+    check_call(('dynamod'+options+' -o '+config).split(), stdout=logfile, stderr=logfile)
+
+    # Run of an equilibration step to blur the system state
+    if ('Rso' in state) and (state['Rso'] != float('inf')) and (state["InitState"] == "Liquid"):
+        print("\n", file=logfile)
+        print("################################", file=logfile)
+        print("#      Liquifaction Run        #", file=logfile)
+        print("################################\n", file=logfile, flush=True)
+        print("# dynarun --unwrapped "+config+" -o "+config+" -c "+str(state['N'] * particle_equil_events)+" --out-data-file data.liqequil.xml.bz2", file=logfile)
+        check_call(["dynarun", "--unwrapped", config, '-o', config, '-c', str(state['N'] * particle_equil_events), "--out-data-file", "data.liqequil.xml.bz2"], stdout=logfile, stderr=logfile)
+
+    # Add the SO Cells global interaction (if needed)
+    if ('Rso' in state) and (state['Rso'] != float('inf')):
+        xml = pydynamo.ConfigFile(config)
+        XMLGlobals = xml.tree.find(".//Globals")
+        XMLSOCells = ET.SubElement(XMLGlobals, 'Global')
+        XMLSOCells.attrib['Name'] = "SOCells" #Name can be anything
+        XMLSOCells.attrib['Type'] = "SOCells" #This must be the right type of Global to load
+        XMLSOCellsRange = ET.SubElement(XMLSOCells, 'Range')
+        XMLSOCellsRange.attrib["Type"] = "All"
+        XMLSOCells.attrib['Diameter'] = str(2 * state['Rso'])
+        xml.save(config)
+
+
+################################################################
+###      DEFINE THE "STATE" VARIABLES TO BE SWEPT & RANGE
+################################################################
+# This is the list of state variables and their ranges
+
+
+statevars = [
+    [ #Sweep 
+        ("Lambda", [float('inf')]),
+        ("InitState", ["FCC"]),
+        ("N", list(map(lambda x: 4*x**3, [3,4,5,6,7]))), #15
+        ('ndensity', list(set(map(lambda x : pydynamo.roundSF(x, 3), [0.01, 0.1, 0.5, 1.0, 1.3])))),
+        ("kT", [1.0]),
+    ],
+]
+
+################################################################
+###          CREATE A SIMULATION MANAGER
+################################################################
+mgr = pydynamo.SimManager("HS_stats", #Which subdirectory to work in
+                          statevars, #State variables
+                          ["p", "CollisionMatrix"], # 'RadialDist' "VACF",  # Output properties
+                          restarts=1, #How many restarts (new initial configurations) should be done per state point
+                          processes=None, #None is automatically use all processors
+)
+
+################################################################
+###          REORGANISE ANY EXISTING SIMULATIONS
+################################################################
+#mgr.reorg_dirs()
+
+################################################################
+###          RUN SOME SIMULATIONS
+################################################################
+mgr.run(setup_worker=setup_worker,
+        particle_equil_events = 1000, # How many events per particle to equilibrate each sim for
+        particle_run_events = 1000, # How many events per particle to run IN TOTAL
+        particle_run_events_block_size=1000) # How big a block each run should be (for jacknife averaging).
+
+################################################################
+###          GET THE DATA
+################################################################
+# This creates a pandas dataframe with columns for the state variables
+# AND any output values. It also generates pkl files, some for
+# different properties.
+df, state_data = mgr.fetch_data(1000)
+
+print(df)

scripts/SW_eos.py

@@ -151,6 +151,6 @@ def setup_worker( config,                # The name of the config file to genera
 # This creates a pandas dataframe with columns for the state variables
 # AND any output values. It also generates pkl files, some for
 # different properties.
-data = mgr.fetch_data(1000)
+df, state_data = mgr.fetch_data(1000)
 
-print(data)
+print(df)