Temperature settings (#98)

* temperature stuff and some input checks

examples/cp2k/H2/infretis.toml

@@ -9,7 +9,7 @@ wmdrun = ['0',
 
 [simulation]
 interfaces = [ 3.45, 3.625, 3.777, 3.922, 4.068, 4.224, 4.399, 4.612, 4.915, 12.0]
-steps = 1000
+steps = 100
 seed = 0
 load_dir = 'load'
 shooting_moves = ['sh','sh', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf']
@@ -21,7 +21,6 @@ allowmaxlength = false
 zero_momentum = true # we do not remove momentum by default in lammps
 n_jumps = 3
 interface_cap = 5.4
-temperature = 300
 
 [engine]
 class = 'cp2k'
@@ -30,6 +29,7 @@ input_path = 'cp2k_input'
 cp2k = 'cp2k.ssmp'
 timestep = 0.2
 subcycles = 10
+temperature = 300
 
 [orderparameter]
 class = 'Distance'

examples/gromacs/H2/gromacs_input/grompp.mdp

@@ -1,4 +1,4 @@
-integrator              = sd
+integrator              = sd ; LGTM
 nsteps                  = 10000
 dt                      = 0.0002
 
@@ -22,13 +22,11 @@ coulombtype             = cutoff
 rcoulomb                = 1.20
 
 tcoupl                  = no
-ref-t			= 300
 tau-t			= 0.1
 tc-grps			= System
 
 pbc                     = xyz
 
-gen_temp                = 300
 gen_seed                = -1
 gen_vel = no
 nstcalcenergy = 1

examples/gromacs/H2/infretis.toml

@@ -1,7 +1,7 @@
 # infretis config
 
 [dask]
-workers = 4
+workers = 1
 wmdrun = ['gmx mdrun -ntomp 1 -ntmpi 1 -pinoffset 0 -pin on',
 	  'gmx mdrun -ntomp 1 -ntmpi 1 -pinoffset 1 -pin on',
 	  'gmx mdrun -ntomp 1 -ntmpi 1 -pinoffset 2 -pin on',
@@ -21,10 +21,6 @@ allowmaxlength = false
 zero_momentum = true # momentum true
 n_jumps = 3
 interface_cap = 0.54
-temperature = 300
-infretis_genvel = true # generate velocities with infretis
-mass = [1.008, 1.008] # then we also need to give masses
-
 
 [engine]
 class = 'gromacs'
@@ -34,6 +30,9 @@ gmx_format = 'g96'
 input_path = 'gromacs_input'
 gmx = 'gmx' # gromacs executable, may also be gmx_mpi, gmx_2023.3, etc.
 subcycles = 10
+temperature = 300
+infretis_genvel = true # generate velocities with infretis
+masses = [1.008, 1.008] # or "masses.txt" stored in input_path/masses.txt in any csv format
 
 [orderparameter]
 class = 'Distance'

examples/gromacs/puckering/README.md

@@ -162,7 +162,7 @@ Fire off `mdrun`. This should take a couple of minutes.
 As you may have guessed by now, a good order parameter for the transition we want to study is the $\theta$ angle. To calculate the angle during the MD run, open `infretis.toml` and replace the indices with the ones you wrote down earlier. You can then recalculate the order parameter using:
 
 ```bash
-inft recalculate_order -trr md.trr -toml infretis.toml -out md-order.txt
+inft recalculate_order -traj md.trr -toml infretis.toml -out md-order.txt
 
 ```
 Plot the $\theta$ values (column 1) vs time (column 0) from the MD run using e.g. gnuplot.

examples/gromacs/puckering/gromacs_input/grompp.mdp

@@ -19,7 +19,6 @@ vdwtype			= cut-off	; Straight cut off for vdw interactions
 tcoupl			= V-rescale	; Canonical sampling thorugh stochastic velocity rescaling
 tc-grps 		= System 	; two coupling groups - more accurate
 tau_t 			= 0.1 		; time constant, in ps
-ref_t 			= 300 		; reference temperature
 
 ; Pressure coupling
 pcoupl 			= no		             ; pressure coupling is on for NPT

examples/gromacs/puckering/step3_infretis/infretis.toml

@@ -24,8 +24,6 @@ shooting_moves = ['sh', 'sh', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf', 'wf', 'w
 maxlength = 2000 # maximum path length (NB! this should never be reached )
 allowmaxlength = false  # do not allow paths to exceed maxlength when shooting from load paths
 zero_momentum = true # remove the center of mass motion when generating vels, done by gromacs by default
-temperature = 300 # velocities are generated at this temperature
-high_accept = true # we do want high-acceptance
 n_jumps = 2 # number of wirefencing jumps
 interface_cap = 70.0 # wirefencing interface cap
 
@@ -37,6 +35,7 @@ gmx_format = 'g96' # only g96 is supported for gromacs
 input_path = '../gromacs_input' # path to topolofy files
 gmx = 'gmx' # gromacs executable, may also be gmx_mpi, gmx_2023.3, etc.
 subcycles = 3 # number of MD steps for each infretis step
+temperature = 300 # temperature of the simulation and velocity generations
 
 [orderparameter]
 class = 'puckering'

examples/lammps/H2/infretis.toml

@@ -1,7 +1,7 @@
 # infretis config
 
 [dask]
-workers = 4
+workers = 1
 wmdrun = ['0',
 	  '0',
 	  '0',
@@ -21,7 +21,6 @@ allowmaxlength = false
 zero_momentum = true # we do not remove momentum by default in lammps
 n_jumps = 3
 interface_cap = 5.4
-temperature = 300
 
 [engine]
 class = 'lammps'
@@ -30,6 +29,7 @@ input_path = 'lammps_input'
 lmp = 'lmp_mpi'
 timestep = 0.2
 subcycles = 10
+temperature = 300
 
 [orderparameter]
 class = 'Distance'

infretis/classes/engines/cp2k.py

@@ -689,6 +689,7 @@ def __init__(
         input_path: str,
         timestep: float,
         subcycles: int,
+        temperature: float,
         extra_files: list[str] | None = None,
         exe_path: str | Path = Path(".").resolve(),
         sleep: float = 0.1,
@@ -700,6 +701,7 @@ def __init__(
             input_path: The path to the directory containing CP2K input files.
             timestep: The time step used in the CP2K simulation.
             subcycles: The number of steps each CP2K run is composed of.
+            temperature: The temperature of the simulation.
             extra_files: List of extra files which may be required to run CP2K.
             exe_path: The path on which the engine is executed
             sleep: A time in seconds, used to wait for files to be ready.
@@ -731,18 +733,27 @@ def __init__(
             ]
         self.mass = np.reshape(mass, (len(mass), 1))
 
-        # read temperature from CP2K input, defaults to 300
-        self.temperature = None
+        # Check that the temperature set in cp2k.inp and infretis.toml
+        # are the same. Often one modifies one and forgets about the other
+        # If no temperature is set in cp2k.inp its NVE with NVT shooting
+        # and we proceed without any remarks (?)
         section = "MOTION->MD"
         nodes = read_cp2k_input(self.input_files["template"])
         node_ref = set_parents(nodes)
         md_settings = node_ref[section]
         for data in md_settings.data:
             if "temperature" in data.lower():
-                self.temperature = float(data.split()[-1])
-        if self.temperature is None:
-            logger.info("No temperature specified in CP2K input. Using 300 K.")
-            self.temperature = 300.0
+                cp2k_temp = float(data.split()[1])
+                if cp2k_temp != temperature:
+                    msg = (
+                        f"The temperature in the cp2k input {cp2k_temp} !="
+                        + f" {temperature} which is specified in the .toml. "
+                        + "They should be the same since we generate "
+                        + "velocities internally at the .toml temperature."
+                    )
+                    raise ValueError(msg)
+
+        self.temperature = temperature
         self.kb = 3.16681534e-6  # hartree
         self.beta = 1 / (self.temperature * self.kb)
 

infretis/classes/engines/enginebase.py

@@ -446,7 +446,7 @@ def _modify_input(
 
     @staticmethod
     def _read_input_settings(
-        sourcefile: str, delim: str = "="
+        sourcefile: str | Path, delim: str = "="
     ) -> dict[str, Any]:
         """Read input settings for simulation input files.
 

infretis/classes/engines/gromacs.py

@@ -92,11 +92,14 @@ def __init__(
         input_path: str | Path,
         timestep: float,
         subcycles: int,
+        temperature: float,
         exe_path: str | Path = Path(".").resolve(),
         maxwarn: int = 0,
         gmx_format: str = "g96",
         write_vel: bool = True,
         write_force: bool = False,
+        infretis_genvel: bool = False,
+        masses: bool | list | str = False,
     ):
         """Set up the GROMACS engine.
 
@@ -107,10 +110,17 @@ def __init__(
             subcycles: The number of steps each GROMACS MD run is composed of.
             exe_path: The absolute path at which the main simulation will be
                 run.
+            temperature: The temperature during the MD run, used when
+                generating velocities and thermostatting by gromacs.
             maxwarn: Setting for the GROMACS `grompp -maxwarn` option.
             gmx_format: The format used for GROMACS configurations.
             write_vel: Determines if GROMACS should write velocities or not.
             write_force: Determines if GROMACS should write forces or not.
+            infretis_genvel: If true we generate the velocities
+                internally by drawing random numbers from a maxwell-boltzmann
+                distribution. Mostly used for testing purposes.
+            masses: A list of particle masses or a txt file with particle
+                masses
         """
         super().__init__("GROMACS engine zamn", timestep, subcycles)
         self.ext = gmx_format
@@ -127,7 +137,7 @@ def __init__(
         # Define the energy terms, these are hard-coded, but
         # here we open up for changing that:
         self.energy_terms = self.select_energy_terms("path")
-        self.input_path = Path(exe_path) / input_path
+        self.input_path = (Path(exe_path) / input_path).resolve()
         # Set the defaults input files:
         default_files = {
             "conf": f"conf.{self.ext}",
@@ -149,10 +159,55 @@ def __init__(
         )
         # Check the input file and create new input file with
         # consistent settings:
+        check_set = self._read_input_settings(self.input_files["input_o"])
+        for key in check_set.keys():
+            val = check_set[key]
+            if (
+                key == "integrator"
+                and "md-vv" not in val
+                and "LGTM" not in val
+            ):
+                msg = (
+                    "Gromacs integrator should be md-vv with "
+                    + "path-sampling. Add a coment LGTM to the relevant"
+                    + " line to overwrite this error."
+                )
+                raise ValueError(msg)
+
+            if key in ["ref-t", "ref_t", "gen-temp", "gen_temp"]:
+                msg = (
+                    f"Found key '{key}' in {self.input_files['input_o']}. "
+                    + "This key is set by infretis and should be removed "
+                    + "from this file."
+                )
+                raise ValueError(msg)
+
+        self.temperature = temperature
+        self.kb = 0.0083144621  # kJ/(K*mol)
+        self.beta = 1 / (self.temperature * self.kb)
+
+        # get masses if generating velocites internally
+        self.infretis_genvel = infretis_genvel
+        if infretis_genvel:
+            if not masses:
+                msg = (
+                    "If you want to generate velocites internally you also "
+                    + "need to specify the particle masses either as a list "
+                    + f"or string pointing to a csv file in {self.input_path}"
+                )
+                raise ValueError(msg)
+            if isinstance(masses, str):
+                self.masses = np.reshape(
+                    np.genfromtxt(self.input_path / masses), (-1, 1)
+                )
+            elif isinstance(masses, list):
+                self.masses = np.reshape(masses, (-1, 1))
+
         settings: dict[str, str | float | int] = {
             "dt": self.timestep,
             "nstxout-compressed": 0,
             "gen_vel": "no",
+            "ref-t": self.temperature,
         }
         for key in (
             "nsteps",
@@ -169,9 +224,6 @@ def __init__(
         if not write_force:
             settings["nstfout"] = 0
 
-        # Add boltzmann constant
-        self.kb = 0.0083144621  # kJ/(K*mol)
-
         # Create the .mdp file name:
         self.input_files["input"] = os.path.join(
             self.input_path, "infretis.mdp"
@@ -639,6 +691,7 @@ def _prepare_shooting_point(
                 "gen_seed": -1,
                 "nsteps": 0,
                 "continuation": "no",
+                "gen-temp": self.temperature,
             }
             self._modify_input(
                 self.input_files["input"], gen_mdp, settings, delim="="
@@ -712,11 +765,6 @@ def modify_velocities(
             vel_settings: A dict containing
                 'zero_momentum': boolean, if true we reset the linear momentum
                   to zero after generating velocities internally.
-                'infretis_genvel': boolen, if true we generate the velocities
-                  internally by drawing random numbers from a maxwell-boltzmann
-                  distribution. Mostly used for testing purposes.
-                'mass': list, the particle masses when generating velocities
-                  internally.
 
 
         Returns:
@@ -728,7 +776,7 @@ def modify_velocities(
         kin_old = system.ekin
         pos = self.dump_frame(system)
         # generate velocities with gromacs
-        if not vel_settings.get("infretis_genvel", False):
+        if not self.infretis_genvel:
             if vel_settings.get("zero_momentum", True) is False:
                 msg = (
                     "Velocitiy generation with gromacs "
@@ -744,18 +792,18 @@ def modify_velocities(
 
         # generate velocities by drawing random numbers
         else:
-            mass = np.reshape(vel_settings["mass"], (-1, 1))
-            beta = 1 / (vel_settings["temperature"] * self.kb)
             txt, xyz, vel, _ = read_gromos96_file(pos)
             if not txt["VELOCITY"]:
                 print(f"{pos} did not contain velocity information.")
                 txt["VELOCITY"] = txt["POSITION"]
-            vel, _ = self.draw_maxwellian_velocities(vel, mass, beta)
+            vel, _ = self.draw_maxwellian_velocities(
+                vel, self.masses, self.beta
+            )
             if vel_settings.get("zero_momentum", False):
-                vel = reset_momentum(vel, mass)
+                vel = reset_momentum(vel, self.masses)
             conf_out = os.path.join(self.exe_dir, f"genvel.{self.ext}")
             write_gromos96_file(conf_out, txt, xyz, vel)
-            kin_new = kinetic_energy(vel, mass)[0]
+            kin_new = kinetic_energy(vel, self.masses)[0]
             system.config = (conf_out, 0)
             system.ekin = kin_new