From 4f310f860884ad7b9e5b1f18e1f3ee7718894956 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Mon, 15 Jul 2024 12:14:33 +0200 Subject: [PATCH 01/23] Start working on Gromacs MD settings --- .../protocols/gromacs_md/md_settings.py | 279 ++++++++++++++++++ 1 file changed, 279 insertions(+) create mode 100644 openfe_gromacs/protocols/gromacs_md/md_settings.py diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py new file mode 100644 index 0000000..2a80375 --- /dev/null +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -0,0 +1,279 @@ +# This code is part of OpenFE and is licensed under the MIT license. +# For details, see https://github.com/OpenFreeEnergy/openfe-gromacs + +"""Settings class for plain MD Protocols using Gromacs + OpenMMTools + +This module implements the settings necessary to run MD simulations using +:class:`openfe.protocols.gromacs_md.md_methods.py` + +""" +from typing import Optional, Literal +from openff.units import unit +from openff.models.types import FloatQuantity +from openfe.protocols.openmm_utils.omm_settings import ( + Settings, + OpenMMSolvationSettings, + OpenMMEngineSettings, + IntegratorSettings, + OpenFFPartialChargeSettings, +) +from gufe.settings import ( + SettingsBaseModel, + OpenMMSystemGeneratorFFSettings +) +try: + from pydantic.v1 import validator +except ImportError: + from pydantic import validator # type: ignore[assignment] + +class MDSimulationSettings(SettingsBaseModel): + """ + Settings for MD simulations in Gromacs. + """ + class Config: + arbitrary_types_allowed = True + + ### Run control ### + integrator: Literal['md', 'md-vv', 'md-vv-avek', 'sd', 'steep'] = 'sd' + """ + MD integrators and other algorithms (steep for energy minimization). + Allowed values are: + 'md': a leap-frog integrator + 'md-vv': a velocity Verlet integrator + 'md-vv-avek': A velocity Verlet algorithm identical to integrator=md-vv, + except that the kinetic energy is determined as the average + of the two half step kinetic energies as in the + integrator=md integrator, and this thus more accurate. + 'sd': A leap-frog stochastic dynamics integrator. Note that when using this + integrator the parameters tcoupl and nsttcouple are ignored. + 'steep': A steepest descent algorithm for energy minimization. + """ + tinit: FloatQuantity['picosecond'] = 0 * unit.picosecond + """ + Starting time for the MD run. + This only makes sense for time-based integrators. + Default 0 * unit.picosecond + """ + dt: FloatQuantity['picosecond'] = 0.001 * unit.picosecond + """ + Time step for integration (only makes sense for time-based integrators). + Default 0.001 * unit.picosecond + """ + nsteps: int = 0 + """ + Maximum number of steps to integrate or minimize, -1 is no maximum + Default 0 + """ + init-step: int = 0 + """ + The starting step. The time at step i in a run is calculated as: + t = tinit + dt * (init-step + i). + Default 0 + """ + simulation-part: int=0 + """ + A simulation can consist of multiple parts, each of which has a part number. + This option specifies what that number will be, which helps keep track of + parts that are logically the same simulation. This option is generally + useful to set only when coping with a crashed simulation where files were + lost. Default 0 + """ + mass-repartition-factor: int = 1 + """ + Scales the masses of the lightest atoms in the system by this factor to + the mass mMin. All atoms with a mass lower than mMin also have their mass + set to that mMin. Default 1 (no mass scaling) + """ + comm-mode: Literal['Linear', 'Angular', 'Linear-acceleration-correction', 'None'] = 'Linear' + """ + Settings for center of mass treatmeant + Allowed values are: + 'Linear': Remove center of mass translational velocity + 'Angular': Remove center of mass translational and rotational velocity + 'Linear-acceleration-correction': Remove center of mass translational + velocity. Correct the center of mass position assuming linear acceleration + over nstcomm steps. + 'None': No restriction on the center of mass motion + """ + nstcomm: int = 100 + """ + Frequency for center of mass motion removal in unit [steps]. + Default 100 + """ + comm-grps: str + """ + Group(s) for center of mass motion removal, default is the whole system. + """ + + ### Langevin dynamics ### + ld-seed: int = -1 + """ + Integer used to initialize random generator for thermal noise for + stochastic and Brownian dynamics. When ld-seed is set to -1, + a pseudo random seed is used. Default -1. + """ + + ### Energy minimization ### + emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) + """ + The minimization is converged when the maximum force is smaller than this + value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) + """ + emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer + """ + Initial step size. Default 0.01 * unit.nanometer + """ + + ### Neighbor searching ### + cutoff-scheme: Literal['verlet'] = 'verlet' + """ + Only allowed option: + 'verlet': Generate a pair list with buffering. The buffer size is + automatically set based on verlet-buffer-tolerance, unless this + is set to -1, in which case rlist will be used. + """ + nstlist: int = 10 + """ + Frequency to update the neighbor list. When dynamics and + verlet-buffer-tolerance set, nstlist is actually a minimum value and + gmx mdrun might increase it, unless it is set to 1. + If set to zero the neighbor list is only constructed once and never updated. + Default 10. + """ + pbc: Literal['xyz', 'no', 'xy'] = 'xyz' + """ + Treatment of periodic boundary conditions. + Allowed values are: + 'xyz': Use periodic boundary conditions in all directions. + 'no': Use no periodic boundary conditions, ignore the box. + 'xy': Use periodic boundary conditions in x and y directions only. + This can be used in combination with walls. + Default 'xyz'. + """ + verlet-buffer-tolerance: FloatQuantity['kilojoule / (mole * picosecond)'] = 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) + """" + Used when performing a simulation with dynamics. This sets the maximum + allowed error for pair interactions per particle caused by the Verlet + buffer, which indirectly sets rlist. + Default 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) + """ + verlet-buffer-pressure-tolerance: FloatQuantity['bar'] = 0.5 * unit.bar + """ + Used when performing a simulation with dynamics and only active when + verlet-buffer-tolerance is positive. This sets the maximum tolerated error + in the average pressure due to missing Lennard-Jones interactions of + particle pairs that are not in the pair list, but come within rvdw range + as the pair list ages. + Default 0.5 * unit.bar + """ + rlist: FloatQuantity['nanometer'] = 1 * unit.nanometer + """ + Cut-off distance for the short-range neighbor list. With dynamics, this is + by default set by the verlet-buffer-tolerance and verlet-buffer-pressure-tolerance + options and the value of rlist is ignored. + """ + + ### Electrostatics ## + + +class MDOutputSettings(SettingsBaseModel): + """" + Output Settings for simulations run using Gromacs + """ + nstxout: int = 0 + """ + Number of steps that elapse between writing coordinates to the output + trajectory file (trr), the last coordinates are always written unless 0, + which means coordinates are not written into the trajectory file. + Default 0. + """ + nstvout: int = 0 + """ + Number of steps that elapse between writing velocities to the output + trajectory file (trr), the last velocities are always written unless 0, + which means velocities are not written into the trajectory file. + Default 0. + """ + nstfout: int = 0 + """ + Number of steps that elapse between writing forces to the output trajectory + file (trr), the last forces are always written, unless 0, which means + forces are not written into the trajectory file. + Default 0. + """ + nstlog: int = 1000 + """ + Number of steps that elapse between writing energies to the log file, the + last energies are always written. Default 1000. + """ + nstcalcenergy: int = 100 + """ + Number of steps that elapse between calculating the energies, 0 is never. + This option is only relevant with dynamics. This option affects the + performance in parallel simulations, because calculating energies requires + global communication between all processes which can become a bottleneck at + high parallelization. Default 100. + """ + nstenergy: int = 1000 + """" + Number of steps that elapse between writing energies to energy file, the + last energies are always written, should be a multiple of nstcalcenergy. + Note that the exact sums and fluctuations over all MD steps modulo + nstcalcenergy are stored in the energy file, so gmx energy can report + exact energy averages and fluctuations also when nstenergy > 1 + """ + nstxout-compressed: int = 0 + """ + Number of steps that elapse between writing position coordinates using + lossy compression (xtc file), 0 for not writing compressed coordinates + output. Default 0. + """ + compressed-x-precision: int = 1000 + """ + Precision with which to write to the compressed trajectory file. + Default 1000. + """ + compressed-x-grps: str + """ + Group(s) to write to the compressed trajectory file, by default the whole + system is written (if nstxout-compressed > 0). + """ + energygrps: str + """ + Group(s) for which to write to write short-ranged non-bonded potential + energies to the energy file (not supported on GPUs) + """ + + +class GromacsMDProtocolSettings(Settings): + class Config: + arbitrary_types_allowed = True + + protocol_repeats: int + """ + Number of independent MD runs to perform. + """ + + @validator('protocol_repeats') + def must_be_positive(cls, v): + if v <= 0: + errmsg = f"protocol_repeats must be a positive value, got {v}." + raise ValueError(errmsg) + return v + + # Things for creating the systems + forcefield_settings: OpenMMSystemGeneratorFFSettings + partial_charge_settings: OpenFFPartialChargeSettings + solvation_settings: OpenMMSolvationSettings + + # MD Engine things + engine_settings: OpenMMEngineSettings + + # Sampling State defining things + integrator_settings: IntegratorSettings + + # Simulation run settings + simulation_settings: MDSimulationSettings + + # Simulations output settings + output_settings: MDOutputSettings From aca3271afa33c52f1fcc89422e6653b9e306a934 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Mon, 15 Jul 2024 15:03:35 +0200 Subject: [PATCH 02/23] adding more settings --- .../protocols/gromacs_md/md_settings.py | 419 +++++++++++++++++- 1 file changed, 418 insertions(+), 1 deletion(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 2a80375..708cb54 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -174,7 +174,424 @@ class Config: """ ### Electrostatics ## - + coulombtype: Literal['cut-off', 'ewald', 'PME', 'P3M-AD', 'Reaction-Field'] = 'PME' + """ + Treatment of electrostatics + Allowed values are: + 'cut-off': Plain cut-off with pair list radius rlist and Coulomb cut-off + rcoulomb, where rlist >= rcoulomb. + 'ewald': Classical Ewald sum electrostatics. + The real-space cut-off rcoulomb should be equal to rlist. + 'PME': Fast smooth Particle-Mesh Ewald (SPME) electrostatics. Direct space + is similar to the Ewald sum, while the reciprocal part is performed + with FFTs. Grid dimensions are controlled with fourierspacing and + the interpolation order with pme-order. + 'P3M-AD': Particle-Particle Particle-Mesh algorithm with analytical + derivative for for long range electrostatic interactions. The + method and code is identical to SPME, except that the influence + function is optimized for the grid. This gives a slight increase + in accuracy. + 'Reaction-Field': Reaction field electrostatics with Coulomb cut-off + rcoulomb, where rlist >= rvdw. The dielectric constant beyond the + cut-off is epsilon-rf. The dielectric constant can be set to + infinity by setting epsilon-rf =0. + Default 'PME' + """ + coulomb-modifier: Literal['Potential-shift', 'None'] = 'Potential-shift' + """ + Allowed options are: + 'Potential-shift': Shift the Coulomb potential by a constant such that it + is zero at the cut-off. + 'None': Use an unmodified Coulomb potential. This can be useful when + comparing energies with those computed with other software. + Default 'Potential-shift' + """ + rcoulomb-switch: FloatQuantity['nanometer'] = 0 * unit.nanometer + """ + Where to start switching the Coulomb potential, only relevant when force + or potential switching is used. + Default 0 * unit.nanometer + """ + rcoulomb: FloatQuantity['nanometer'] = 1.2 * unit.nanometer + """" + The distance for the Coulomb cut-off. Note that with PME this value can be + increased by the PME tuning in gmx mdrun along with the PME grid spacing. + Default 1.2 * unit.nanometer + """ + epsilon-r: int = 1 + """ + The relative dielectric constant. A value of 0 means infinity. Default 1. + """ + epsilon-rf: int = 0 + """ + The relative dielectric constant of the reaction field. This is only used + with reaction-field electrostatics. A value of 0 means infinity. Default 0 + """ + + ### Van der Waals ### + vdwtype: Literal['Cut-off', 'PME'] = 'Cut-off' + """ + Treatment of vdW interactions. Allowed options are: + 'Cut-off': Plain cut-off with pair list radius rlist and VdW cut-off rvdw, + where rlist >= rvdw. + 'PME: Fast smooth Particle-mesh Ewald (SPME) for VdW interactions. + Default 'Cut-off'. + """ + vdw-modifier: Literal['Potential-shift', 'None', 'Force-switch', 'Potential-switch'] = 'Potential-shift' + """ + Allowed values are: + 'Potential-shift': Shift the Van der Waals potential by a constant such + that it is zero at the cut-off. + 'None': Use an unmodified Van der Waals potential. This can be useful when + comparing energies with those computed with other software. + 'Force-switch': Smoothly switches the forces to zero between rvdw-switch + and rvdw. This shifts the potential shift over the whole range and + switches it to zero at the cut-off. Note that this is more expensive + to calculate than a plain cut-off and it is not required for energy + conservation, since Potential-shift conserves energy just as well. + 'Potential-switch': Smoothly switches the potential to zero between + rvdw-switch and rvdw. Note that this introduces articifically large + forces in the switching region and is much more expensive to calculate. + This option should only be used if the force field you are using + requires this. + """ + rvdw-switch: FloatQuantity['nanometer'] = 0 * unit.nanometer + """" + Where to start switching the LJ force and possibly the potential, only + relevant when force or potential switching is used. + Default 0 * unit.nanometer + """ + rvdw: FloatQuantity['nanometer'] = 1.0 * unit.nanometer + """ + Distance for the LJ or Buckingham cut-off. Default 1 * unit.nanometer + """ + DispCorr: Literal['no', 'EnerPres', 'Ener'] = 'EnerPres' + """ + Allowed values are: + 'no': Don’t apply any correction + 'EnerPres': Apply long range dispersion corrections for Energy and Pressure. + 'Ener': Apply long range dispersion corrections for Energy only. + Default 'EnerPres' + """ + + ### Ewald ### + fourierspacing: FloatQuantity['nanometer'] = 0.12 * unit.nanometer + """ + For ordinary Ewald, the ratio of the box dimensions and the spacing + determines a lower bound for the number of wave vectors to use in each + (signed) direction. For PME and P3M, that ratio determines a lower bound + for the number of Fourier-space grid points that will be used along that + axis. In all cases, the number for each direction can be overridden by + entering a non-zero value for that fourier-nx direction. + Default 0.12 * unit.nanometer + """ + pme-order: int = 4 + """ + The number of grid points along a dimension to which a charge is mapped. + The actual order of the PME interpolation is one less, e.g. the default of + 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for P3M-AD). + When running in parallel, it can be worth to switch to 5 and simultaneously + increase the grid spacing. Note that on the CPU only values 4 and 5 have + SIMD acceleration and GPUs only support the value 4. + Default 4. + """ + ewald-rtol: float = 1e-5 + """ + The relative strength of the Ewald-shifted direct potential at rcoulomb is + given by ewald-rtol. Decreasing this will give a more accurate direct sum, + but then you need more wave vectors for the reciprocal sum. + Default 1e-5 + """ + lj-pme-comb-rule: Literal['Geometric', 'Lorentz-Berthelot'] = 'Geometric' + """ + The combination rules used to combine VdW-parameters in the reciprocal part + of LJ-PME. Geometric rules are much faster than Lorentz-Berthelot and + usually the recommended choice, even when the rest of the force field uses + the Lorentz-Berthelot rules. Allowed values are: + 'Geometric': Apply geometric combination rules. + 'Lorentz-Berthelot': Apply Lorentz-Berthelot combination rules. + Default 'Geometric'. + """ + ewald-geometry: Literal['3d', '3dc'] = '3d' + """ + Allowed values are: + '3d': The Ewald sum is performed in all three dimensions. + '3dc': The reciprocal sum is still performed in 3D, but a force and + potential correction applied in the z dimension to produce a pseudo-2D summation. + Default '3d'. + """ + epsilon-surface: float = 0 + """ + This controls the dipole correction to the Ewald summation in 3D. + The default value of zero means it is turned off. Turn it on by setting it + to the value of the relative permittivity of the imaginary surface around + your infinite system. Be careful - you shouldn’t use this if you have free + mobile charges in your system. This value does not affect the slab 3DC + variant of the long range corrections. Default 0. + """ + + ### Temerature coupling ### + tcoupl: Literal['no', 'berendsen', 'nose-hoover', 'andersen', 'andersen-massive', 'v-rescale'] = 'no' + """ + Temperature coupling options. Note that tcoupl will be ignored when the + 'sd' integrator is used. + Allowed values are: + 'no': No temperature coupling. + 'berendsen': Temperature coupling with a Berendsen thermostat to a bath + with temperature ref-t, with time constant tau-t. Several groups can be + coupled separately, these are specified in the tc-grps field separated + by spaces. This is a historical thermostat needed to be able to + reproduce previous simulations, but we strongly recommend not to use + it for new production runs. + 'nose-hoover': Temperature coupling using a Nose-Hoover extended ensemble. + The reference temperature and coupling groups are selected as above, + but in this case tau-t controls the period of the temperature + fluctuations at equilibrium, which is slightly different from a + relaxation time. For NVT simulations the conserved energy quantity is + written to the energy and log files. + 'andersen': Temperature coupling by randomizing a fraction of the particle + velocities at each timestep. Reference temperature and coupling groups + are selected as above. tau-t is the average time between randomization + of each molecule. Inhibits particle dynamics somewhat, but little or no + ergodicity issues. Currently only implemented with velocity Verlet, and + not implemented with constraints. + 'andersen-massive': Temperature coupling by randomizing velocities of all + particles at infrequent timesteps. Reference temperature and coupling + groups are selected as above. tau-t is the time between randomization + of all molecules. Inhibits particle dynamics somewhat, but little or + no ergodicity issues. Currently only implemented with velocity Verlet. + 'v-rescale': Temperature coupling using velocity rescaling with a + stochastic term (JCP 126, 014101). This thermostat is similar to + Berendsen coupling, with the same scaling using tau-t, but the + stochastic term ensures that a proper canonical ensemble is generated. + The random seed is set with ld-seed. This thermostat works correctly + even for tau-t =0. For NVT simulations the conserved energy quantity + is written to the energy and log file. + Default 'no' (since for the default integrator, 'sd', this option is ignored). + """ + nsttcouple: int = -1 + """ + The frequency for coupling the temperature. The default value of -1 sets + nsttcouple equal to 100, or fewer steps if required for accurate + integration (5 steps per tau for first order coupling, 20 steps per tau for + second order coupling). Note that the default value is large in order to + reduce the overhead of the additional computation and communication + required for obtaining the kinetic energy. For velocity Verlet integrators + nsttcouple is set to 1. + Default -1. + """ + tc-grps: str = 'system' + """ + Groups to couple to separate temperature baths. Default 'system'. + """ + tau-t: FloatQuantity['picosecond'] = 2.0 * unit.picosecond + """" + Time constant for coupling (one for each group in tc-grps), -1 means no + temperature coupling. Default 2.0 * unit.picosecond. + """ + ref-t: FloatQuantity['Kelvin'] = 298.15 * unit.kelvin + """ + Reference temperature for coupling (one for each group in tc-grps). + Default 298.15 * unit.kelvin + """ + + ### Pressure coupling ### + pcoupl: Literal['no', 'berendsen', 'C-rescale', 'Parrinello-Rahman', 'MTTK'] = 'no' + """ + Options for pressure coupling (barostat). Allowed values are: + 'no': No pressure coupling. This means a fixed box size. + 'berendsen': Exponential relaxation pressure coupling with time constant + tau-p. The box is scaled every nstpcouple steps. This barostat does not + yield a correct thermodynamic ensemble; it is only included to be able + to reproduce previous runs, and we strongly recommend against using it + for new simulations. + 'C-rescale': Exponential relaxation pressure coupling with time constant + tau-p, including a stochastic term to enforce correct volume + fluctuations. The box is scaled every nstpcouple steps. It can be used + for both equilibration and production. + 'Parrinello-Rahman': Extended-ensemble pressure coupling where the box + vectors are subject to an equation of motion. The equation of motion + for the atoms is coupled to this. No instantaneous scaling takes place. + As for Nose-Hoover temperature coupling the time constant tau-p is the + period of pressure fluctuations at equilibrium. + 'MTTK': Martyna-Tuckerman-Tobias-Klein implementation, only useable with + integrator=md-vv or integrator=md-vv-avek, very similar to Parrinello-Rahman. + As for Nose-Hoover temperature coupling the time constant tau-p is the + period of pressure fluctuations at equilibrium. This is probably a + better method when you want to apply pressure scaling during data + collection, but beware that you can get very large oscillations if you + are starting from a different pressure. This requires a constant ensemble + temperature for the system. Currently it only supports isotropic scaling, + and only works without constraints. + Default 'no'. + """ + pcoupltype: Literal['isotropic', 'semiisotropic', 'anisotropic', 'surface-tension'] = 'isotropic' + """ + Specifies the kind of isotropy of the pressure coupling used. Each kind + takes one or more values for compressibility and ref-p. Only a single value + is permitted for tau-p. Allowed values are: + 'isotropic': sotropic pressure coupling with time constant tau-p. One value + each for compressibility and ref-p is required. + 'semiisotropic': Pressure coupling which is isotropic in the x and y + direction, but different in the z direction. This can be useful for + membrane simulations. Two values each for compressibility and ref-p are + required, for x/y and z directions respectively. + 'anisotropic': Same as before, but 6 values are needed for xx, yy, zz, + xy/yx, xz/zx and yz/zy components, respectively. When the off-diagonal + compressibilities are set to zero, a rectangular box will stay + rectangular. Beware that anisotropic scaling can lead to extreme + deformation of the simulation box. + 'surface-tension': Surface tension coupling for surfaces parallel to the + xy-plane. Uses normal pressure coupling for the z-direction, while the + surface tension is coupled to the x/y dimensions of the box. The first + ref-p value is the reference surface tension times the number of + surfaces bar nm, the second value is the reference z-pressure bar. + The two compressibility values are the compressibility in the x/y and z + direction respectively. The value for the z-compressibility should be + reasonably accurate since it influences the convergence of the + surface-tension, it can also be set to zero to have a box with constant + height. + Default 'isotropic' + """ + nstpcouple: int = -1 + """ + The frequency for coupling the pressure. The default value of -1 sets + nstpcouple equal to 100, or fewer steps if required for accurate integration + (5 steps per tau for first order coupling, 20 steps per tau for second order + coupling). Note that the default value is large in order to reduce the + overhead of the additional computation and communication required for + obtaining the virial and kinetic energy. For velocity Verlet integrators + nsttcouple is set to 1. Default -1. + """ + tau-p: FloatQuantity['picosecond'] = 5 * unit.picosecond + """ + The time constant for pressure coupling (one value for all directions). + Default 5 * unit.picosecond. + """ + compressibility: FloatQuantity['1/bar'] = 4.5e-05 / unit.bar + """ + The compressibility. For water at 1 atm and 300 K the compressibility is + 4.5e-5 bar-1. The number of required values is implied by pcoupltype. + Default 4.5e-05 / unit.bar + """ + ref-p: FloatQuantity['bar'] = 1.01325 * unit.bar + """ + The reference pressure for coupling. The number of required values is + implied by pcoupltype. Default 1.01325 * unit.bar. + """ + refcoord-scaling: Literal['no', 'all', 'com'] = 'no' + """ + Allowed values are: + 'no': The reference coordinates for position restraints are not modified. + 'all': The reference coordinates are scaled with the scaling matrix of the + pressure coupling. + 'com': Scale the center of mass of the reference coordinates with the + scaling matrix of the pressure coupling. The vectors of each reference + coordinate to the center of mass are not scaled. Only one COM is used, + even when there are multiple molecules with position restraints. + For calculating the COM of the reference coordinates in the starting + configuration, periodic boundary conditions are not taken into account. + Default 'no'. + """ + + ### Velocity generation ### + gen-vel: Literal['no', 'yes'] = 'yes' + """ + Velocity generation. Allowed values are: + 'no': Do not generate velocities. The velocities are set to zero when there + are no velocities in the input structure file. + 'yes': Generate velocities in gmx grompp according to a Maxwell distribution + at temperature gen-temp, with random seed gen-seed. This is only + meaningful with integrator=md. + Default 'yes'. + """ + gen-temp: FloatQuantity['kelvin'] = 298.15 * unit.kelvin + """ + Temperature for Maxwell distribution. Default 298.15 * unit.kelvin. + """ + gen-seed: int = -1 + """ + Used to initialize random generator for random velocities, when gen-seed is + set to -1, a pseudo random seed is used. Default -1. + """ + + ### Bonds ### + constraints: Literal['none', 'h-bonds', 'all-bonds', 'h-angles', 'all-angles'] = 'h-bonds' + """ + Controls which bonds in the topology will be converted to rigid holonomic + constraints. Note that typical rigid water models do not have bonds, but + rather a specialized [settles] directive, so are not affected by this keyword. + Allowed values are: + 'none': No bonds converted to constraints. + 'h-bonds': Convert the bonds with H-atoms to constraints. + 'all-bonds': Convert all bonds to constraints. + 'h-angles': Convert all bonds to constraints and convert the angles that + involve H-atoms to bond-constraints. + 'all-angles': Convert all bonds to constraints and all angles to + bond-constraints. + Default 'h-bonds' + """ + constraint-algorithm: Literal['lincs', 'shake'] = 'lincs' + """ + Chooses which solver satisfies any non-SETTLE holonomic constraints. + Allowed values are: + 'lincs': LINear Constraint Solver. With domain decomposition the parallel + version P-LINCS is used. The accuracy in set with lincs-order, which + sets the number of matrices in the expansion for the matrix inversion. + After the matrix inversion correction the algorithm does an iterative + correction to compensate for lengthening due to rotation. The number of + such iterations can be controlled with lincs-iter. The root mean square + relative constraint deviation is printed to the log file every nstlog + steps. If a bond rotates more than lincs-warnangle in one step, a + warning will be printed both to the log file and to stderr. + LINCS should not be used with coupled angle constraints. + 'shake': SHAKE is slightly slower and less stable than LINCS, but does work + with angle constraints. The relative tolerance is set with shake-tol, + 0.0001 is a good value for “normal” MD. SHAKE does not support + constraints between atoms on different decomposition domains, so it can + only be used with domain decomposition when so-called update-groups are + used, which is usally the case when only bonds involving hydrogens are + constrained. SHAKE can not be used with energy minimization. + Default 'lincs' + """ + continuation: Literal['no', 'yes'] = 'no' + """ + This option was formerly known as unconstrained-start. + Allowed values are: + 'no': Apply constraints to the start configuration and reset shells. + 'yes': Do not apply constraints to the start configuration and do not reset + shells, useful for exact coninuation and reruns. + """ + shake-tol: float = 0.0001 + """ Relative tolerance for SHAKE. Default 0.0001. """ + lincs-order: int = 12 + """ + Highest order in the expansion of the constraint coupling matrix. When + constraints form triangles, an additional expansion of the same order is + applied on top of the normal expansion only for the couplings within such + triangles. For “normal” MD simulations an order of 4 usually suffices, 6 is + needed for large time-steps with virtual sites or BD. For accurate energy + minimization in double precision an order of 8 or more might be required. + Note that in single precision an order higher than 6 will often lead to + worse accuracy due to amplification of rounding errors. Default 12. + """ + lincs-iter: int = 1 + """ + Number of iterations to correct for rotational lengthening in LINCS. + Default 1. + """ + lincs-warnangle: FloatQuantity['deg'] = 30 * unit.degree + """ + Maximum angle that a bond can rotate before LINCS will complain. + Default 30 * unit.degree. + """ + morse: Literal['no', 'yes'] = 'no' + """ + Allowed options are: + 'no': Bonds are represented by a harmonic potential. + 'yes': Bonds are represented by a Morse potential. + Default 'no'. + """ class MDOutputSettings(SettingsBaseModel): """" From 10d6f0963fe5c092695dbee628812d518a30d37a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 15 Jul 2024 13:04:18 +0000 Subject: [PATCH 03/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../protocols/gromacs_md/md_settings.py | 461 +++++++++--------- 1 file changed, 230 insertions(+), 231 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 708cb54..9b87ead 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -7,20 +7,19 @@ :class:`openfe.protocols.gromacs_md.md_methods.py` """ -from typing import Optional, Literal -from openff.units import unit -from openff.models.types import FloatQuantity +from typing import Literal, Optional + +from gufe.settings import OpenMMSystemGeneratorFFSettings, SettingsBaseModel from openfe.protocols.openmm_utils.omm_settings import ( - Settings, - OpenMMSolvationSettings, - OpenMMEngineSettings, IntegratorSettings, OpenFFPartialChargeSettings, + OpenMMEngineSettings, + OpenMMSolvationSettings, + Settings, ) -from gufe.settings import ( - SettingsBaseModel, - OpenMMSystemGeneratorFFSettings -) +from openff.models.types import FloatQuantity +from openff.units import unit + try: from pydantic.v1 import validator except ImportError: @@ -40,19 +39,19 @@ class Config: Allowed values are: 'md': a leap-frog integrator 'md-vv': a velocity Verlet integrator - 'md-vv-avek': A velocity Verlet algorithm identical to integrator=md-vv, - except that the kinetic energy is determined as the average - of the two half step kinetic energies as in the + 'md-vv-avek': A velocity Verlet algorithm identical to integrator=md-vv, + except that the kinetic energy is determined as the average + of the two half step kinetic energies as in the integrator=md integrator, and this thus more accurate. 'sd': A leap-frog stochastic dynamics integrator. Note that when using this - integrator the parameters tcoupl and nsttcouple are ignored. + integrator the parameters tcoupl and nsttcouple are ignored. 'steep': A steepest descent algorithm for energy minimization. """ tinit: FloatQuantity['picosecond'] = 0 * unit.picosecond """ - Starting time for the MD run. + Starting time for the MD run. This only makes sense for time-based integrators. - Default 0 * unit.picosecond + Default 0 * unit.picosecond """ dt: FloatQuantity['picosecond'] = 0.001 * unit.picosecond """ @@ -66,22 +65,22 @@ class Config: """ init-step: int = 0 """ - The starting step. The time at step i in a run is calculated as: + The starting step. The time at step i in a run is calculated as: t = tinit + dt * (init-step + i). Default 0 """ simulation-part: int=0 """ A simulation can consist of multiple parts, each of which has a part number. - This option specifies what that number will be, which helps keep track of - parts that are logically the same simulation. This option is generally - useful to set only when coping with a crashed simulation where files were + This option specifies what that number will be, which helps keep track of + parts that are logically the same simulation. This option is generally + useful to set only when coping with a crashed simulation where files were lost. Default 0 """ mass-repartition-factor: int = 1 """ - Scales the masses of the lightest atoms in the system by this factor to - the mass mMin. All atoms with a mass lower than mMin also have their mass + Scales the masses of the lightest atoms in the system by this factor to + the mass mMin. All atoms with a mass lower than mMin also have their mass set to that mMin. Default 1 (no mass scaling) """ comm-mode: Literal['Linear', 'Angular', 'Linear-acceleration-correction', 'None'] = 'Linear' @@ -90,14 +89,14 @@ class Config: Allowed values are: 'Linear': Remove center of mass translational velocity 'Angular': Remove center of mass translational and rotational velocity - 'Linear-acceleration-correction': Remove center of mass translational - velocity. Correct the center of mass position assuming linear acceleration - over nstcomm steps. + 'Linear-acceleration-correction': Remove center of mass translational + velocity. Correct the center of mass position assuming linear acceleration + over nstcomm steps. 'None': No restriction on the center of mass motion """ nstcomm: int = 100 """ - Frequency for center of mass motion removal in unit [steps]. + Frequency for center of mass motion removal in unit [steps]. Default 100 """ comm-grps: str @@ -108,15 +107,15 @@ class Config: ### Langevin dynamics ### ld-seed: int = -1 """ - Integer used to initialize random generator for thermal noise for - stochastic and Brownian dynamics. When ld-seed is set to -1, + Integer used to initialize random generator for thermal noise for + stochastic and Brownian dynamics. When ld-seed is set to -1, a pseudo random seed is used. Default -1. """ ### Energy minimization ### emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) """ - The minimization is converged when the maximum force is smaller than this + The minimization is converged when the maximum force is smaller than this value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) """ emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer @@ -128,15 +127,15 @@ class Config: cutoff-scheme: Literal['verlet'] = 'verlet' """ Only allowed option: - 'verlet': Generate a pair list with buffering. The buffer size is - automatically set based on verlet-buffer-tolerance, unless this + 'verlet': Generate a pair list with buffering. The buffer size is + automatically set based on verlet-buffer-tolerance, unless this is set to -1, in which case rlist will be used. """ nstlist: int = 10 """ - Frequency to update the neighbor list. When dynamics and - verlet-buffer-tolerance set, nstlist is actually a minimum value and - gmx mdrun might increase it, unless it is set to 1. + Frequency to update the neighbor list. When dynamics and + verlet-buffer-tolerance set, nstlist is actually a minimum value and + gmx mdrun might increase it, unless it is set to 1. If set to zero the neighbor list is only constructed once and never updated. Default 10. """ @@ -146,30 +145,30 @@ class Config: Allowed values are: 'xyz': Use periodic boundary conditions in all directions. 'no': Use no periodic boundary conditions, ignore the box. - 'xy': Use periodic boundary conditions in x and y directions only. + 'xy': Use periodic boundary conditions in x and y directions only. This can be used in combination with walls. Default 'xyz'. """ verlet-buffer-tolerance: FloatQuantity['kilojoule / (mole * picosecond)'] = 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) """" - Used when performing a simulation with dynamics. This sets the maximum - allowed error for pair interactions per particle caused by the Verlet + Used when performing a simulation with dynamics. This sets the maximum + allowed error for pair interactions per particle caused by the Verlet buffer, which indirectly sets rlist. Default 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) """ verlet-buffer-pressure-tolerance: FloatQuantity['bar'] = 0.5 * unit.bar """ - Used when performing a simulation with dynamics and only active when - verlet-buffer-tolerance is positive. This sets the maximum tolerated error - in the average pressure due to missing Lennard-Jones interactions of - particle pairs that are not in the pair list, but come within rvdw range - as the pair list ages. + Used when performing a simulation with dynamics and only active when + verlet-buffer-tolerance is positive. This sets the maximum tolerated error + in the average pressure due to missing Lennard-Jones interactions of + particle pairs that are not in the pair list, but come within rvdw range + as the pair list ages. Default 0.5 * unit.bar """ rlist: FloatQuantity['nanometer'] = 1 * unit.nanometer """ - Cut-off distance for the short-range neighbor list. With dynamics, this is - by default set by the verlet-buffer-tolerance and verlet-buffer-pressure-tolerance + Cut-off distance for the short-range neighbor list. With dynamics, this is + by default set by the verlet-buffer-tolerance and verlet-buffer-pressure-tolerance options and the value of rlist is ignored. """ @@ -178,43 +177,43 @@ class Config: """ Treatment of electrostatics Allowed values are: - 'cut-off': Plain cut-off with pair list radius rlist and Coulomb cut-off + 'cut-off': Plain cut-off with pair list radius rlist and Coulomb cut-off rcoulomb, where rlist >= rcoulomb. - 'ewald': Classical Ewald sum electrostatics. + 'ewald': Classical Ewald sum electrostatics. The real-space cut-off rcoulomb should be equal to rlist. - 'PME': Fast smooth Particle-Mesh Ewald (SPME) electrostatics. Direct space - is similar to the Ewald sum, while the reciprocal part is performed - with FFTs. Grid dimensions are controlled with fourierspacing and + 'PME': Fast smooth Particle-Mesh Ewald (SPME) electrostatics. Direct space + is similar to the Ewald sum, while the reciprocal part is performed + with FFTs. Grid dimensions are controlled with fourierspacing and the interpolation order with pme-order. - 'P3M-AD': Particle-Particle Particle-Mesh algorithm with analytical - derivative for for long range electrostatic interactions. The - method and code is identical to SPME, except that the influence - function is optimized for the grid. This gives a slight increase + 'P3M-AD': Particle-Particle Particle-Mesh algorithm with analytical + derivative for for long range electrostatic interactions. The + method and code is identical to SPME, except that the influence + function is optimized for the grid. This gives a slight increase in accuracy. - 'Reaction-Field': Reaction field electrostatics with Coulomb cut-off - rcoulomb, where rlist >= rvdw. The dielectric constant beyond the - cut-off is epsilon-rf. The dielectric constant can be set to + 'Reaction-Field': Reaction field electrostatics with Coulomb cut-off + rcoulomb, where rlist >= rvdw. The dielectric constant beyond the + cut-off is epsilon-rf. The dielectric constant can be set to infinity by setting epsilon-rf =0. Default 'PME' """ coulomb-modifier: Literal['Potential-shift', 'None'] = 'Potential-shift' """ Allowed options are: - 'Potential-shift': Shift the Coulomb potential by a constant such that it - is zero at the cut-off. - 'None': Use an unmodified Coulomb potential. This can be useful when + 'Potential-shift': Shift the Coulomb potential by a constant such that it + is zero at the cut-off. + 'None': Use an unmodified Coulomb potential. This can be useful when comparing energies with those computed with other software. Default 'Potential-shift' """ rcoulomb-switch: FloatQuantity['nanometer'] = 0 * unit.nanometer """ - Where to start switching the Coulomb potential, only relevant when force + Where to start switching the Coulomb potential, only relevant when force or potential switching is used. Default 0 * unit.nanometer """ rcoulomb: FloatQuantity['nanometer'] = 1.2 * unit.nanometer """" - The distance for the Coulomb cut-off. Note that with PME this value can be + The distance for the Coulomb cut-off. Note that with PME this value can be increased by the PME tuning in gmx mdrun along with the PME grid spacing. Default 1.2 * unit.nanometer """ @@ -224,7 +223,7 @@ class Config: """ epsilon-rf: int = 0 """ - The relative dielectric constant of the reaction field. This is only used + The relative dielectric constant of the reaction field. This is only used with reaction-field electrostatics. A value of 0 means infinity. Default 0 """ @@ -232,7 +231,7 @@ class Config: vdwtype: Literal['Cut-off', 'PME'] = 'Cut-off' """ Treatment of vdW interactions. Allowed options are: - 'Cut-off': Plain cut-off with pair list radius rlist and VdW cut-off rvdw, + 'Cut-off': Plain cut-off with pair list radius rlist and VdW cut-off rvdw, where rlist >= rvdw. 'PME: Fast smooth Particle-mesh Ewald (SPME) for VdW interactions. Default 'Cut-off'. @@ -240,25 +239,25 @@ class Config: vdw-modifier: Literal['Potential-shift', 'None', 'Force-switch', 'Potential-switch'] = 'Potential-shift' """ Allowed values are: - 'Potential-shift': Shift the Van der Waals potential by a constant such + 'Potential-shift': Shift the Van der Waals potential by a constant such that it is zero at the cut-off. - 'None': Use an unmodified Van der Waals potential. This can be useful when + 'None': Use an unmodified Van der Waals potential. This can be useful when comparing energies with those computed with other software. - 'Force-switch': Smoothly switches the forces to zero between rvdw-switch - and rvdw. This shifts the potential shift over the whole range and - switches it to zero at the cut-off. Note that this is more expensive - to calculate than a plain cut-off and it is not required for energy + 'Force-switch': Smoothly switches the forces to zero between rvdw-switch + and rvdw. This shifts the potential shift over the whole range and + switches it to zero at the cut-off. Note that this is more expensive + to calculate than a plain cut-off and it is not required for energy conservation, since Potential-shift conserves energy just as well. - 'Potential-switch': Smoothly switches the potential to zero between - rvdw-switch and rvdw. Note that this introduces articifically large - forces in the switching region and is much more expensive to calculate. - This option should only be used if the force field you are using + 'Potential-switch': Smoothly switches the potential to zero between + rvdw-switch and rvdw. Note that this introduces articifically large + forces in the switching region and is much more expensive to calculate. + This option should only be used if the force field you are using requires this. """ rvdw-switch: FloatQuantity['nanometer'] = 0 * unit.nanometer """" - Where to start switching the LJ force and possibly the potential, only - relevant when force or potential switching is used. + Where to start switching the LJ force and possibly the potential, only + relevant when force or potential switching is used. Default 0 * unit.nanometer """ rvdw: FloatQuantity['nanometer'] = 1.0 * unit.nanometer @@ -277,36 +276,36 @@ class Config: ### Ewald ### fourierspacing: FloatQuantity['nanometer'] = 0.12 * unit.nanometer """ - For ordinary Ewald, the ratio of the box dimensions and the spacing - determines a lower bound for the number of wave vectors to use in each - (signed) direction. For PME and P3M, that ratio determines a lower bound - for the number of Fourier-space grid points that will be used along that - axis. In all cases, the number for each direction can be overridden by + For ordinary Ewald, the ratio of the box dimensions and the spacing + determines a lower bound for the number of wave vectors to use in each + (signed) direction. For PME and P3M, that ratio determines a lower bound + for the number of Fourier-space grid points that will be used along that + axis. In all cases, the number for each direction can be overridden by entering a non-zero value for that fourier-nx direction. Default 0.12 * unit.nanometer """ pme-order: int = 4 """ - The number of grid points along a dimension to which a charge is mapped. - The actual order of the PME interpolation is one less, e.g. the default of - 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for P3M-AD). - When running in parallel, it can be worth to switch to 5 and simultaneously - increase the grid spacing. Note that on the CPU only values 4 and 5 have - SIMD acceleration and GPUs only support the value 4. + The number of grid points along a dimension to which a charge is mapped. + The actual order of the PME interpolation is one less, e.g. the default of + 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for P3M-AD). + When running in parallel, it can be worth to switch to 5 and simultaneously + increase the grid spacing. Note that on the CPU only values 4 and 5 have + SIMD acceleration and GPUs only support the value 4. Default 4. """ ewald-rtol: float = 1e-5 """ - The relative strength of the Ewald-shifted direct potential at rcoulomb is - given by ewald-rtol. Decreasing this will give a more accurate direct sum, + The relative strength of the Ewald-shifted direct potential at rcoulomb is + given by ewald-rtol. Decreasing this will give a more accurate direct sum, but then you need more wave vectors for the reciprocal sum. Default 1e-5 """ lj-pme-comb-rule: Literal['Geometric', 'Lorentz-Berthelot'] = 'Geometric' """ - The combination rules used to combine VdW-parameters in the reciprocal part - of LJ-PME. Geometric rules are much faster than Lorentz-Berthelot and - usually the recommended choice, even when the rest of the force field uses + The combination rules used to combine VdW-parameters in the reciprocal part + of LJ-PME. Geometric rules are much faster than Lorentz-Berthelot and + usually the recommended choice, even when the rest of the force field uses the Lorentz-Berthelot rules. Allowed values are: 'Geometric': Apply geometric combination rules. 'Lorentz-Berthelot': Apply Lorentz-Berthelot combination rules. @@ -316,67 +315,67 @@ class Config: """ Allowed values are: '3d': The Ewald sum is performed in all three dimensions. - '3dc': The reciprocal sum is still performed in 3D, but a force and + '3dc': The reciprocal sum is still performed in 3D, but a force and potential correction applied in the z dimension to produce a pseudo-2D summation. Default '3d'. """ epsilon-surface: float = 0 """ - This controls the dipole correction to the Ewald summation in 3D. - The default value of zero means it is turned off. Turn it on by setting it - to the value of the relative permittivity of the imaginary surface around - your infinite system. Be careful - you shouldn’t use this if you have free - mobile charges in your system. This value does not affect the slab 3DC + This controls the dipole correction to the Ewald summation in 3D. + The default value of zero means it is turned off. Turn it on by setting it + to the value of the relative permittivity of the imaginary surface around + your infinite system. Be careful - you shouldn’t use this if you have free + mobile charges in your system. This value does not affect the slab 3DC variant of the long range corrections. Default 0. """ ### Temerature coupling ### tcoupl: Literal['no', 'berendsen', 'nose-hoover', 'andersen', 'andersen-massive', 'v-rescale'] = 'no' """ - Temperature coupling options. Note that tcoupl will be ignored when the + Temperature coupling options. Note that tcoupl will be ignored when the 'sd' integrator is used. - Allowed values are: + Allowed values are: 'no': No temperature coupling. - 'berendsen': Temperature coupling with a Berendsen thermostat to a bath - with temperature ref-t, with time constant tau-t. Several groups can be - coupled separately, these are specified in the tc-grps field separated - by spaces. This is a historical thermostat needed to be able to - reproduce previous simulations, but we strongly recommend not to use + 'berendsen': Temperature coupling with a Berendsen thermostat to a bath + with temperature ref-t, with time constant tau-t. Several groups can be + coupled separately, these are specified in the tc-grps field separated + by spaces. This is a historical thermostat needed to be able to + reproduce previous simulations, but we strongly recommend not to use it for new production runs. - 'nose-hoover': Temperature coupling using a Nose-Hoover extended ensemble. - The reference temperature and coupling groups are selected as above, - but in this case tau-t controls the period of the temperature - fluctuations at equilibrium, which is slightly different from a - relaxation time. For NVT simulations the conserved energy quantity is + 'nose-hoover': Temperature coupling using a Nose-Hoover extended ensemble. + The reference temperature and coupling groups are selected as above, + but in this case tau-t controls the period of the temperature + fluctuations at equilibrium, which is slightly different from a + relaxation time. For NVT simulations the conserved energy quantity is written to the energy and log files. - 'andersen': Temperature coupling by randomizing a fraction of the particle - velocities at each timestep. Reference temperature and coupling groups - are selected as above. tau-t is the average time between randomization - of each molecule. Inhibits particle dynamics somewhat, but little or no - ergodicity issues. Currently only implemented with velocity Verlet, and + 'andersen': Temperature coupling by randomizing a fraction of the particle + velocities at each timestep. Reference temperature and coupling groups + are selected as above. tau-t is the average time between randomization + of each molecule. Inhibits particle dynamics somewhat, but little or no + ergodicity issues. Currently only implemented with velocity Verlet, and not implemented with constraints. - 'andersen-massive': Temperature coupling by randomizing velocities of all - particles at infrequent timesteps. Reference temperature and coupling - groups are selected as above. tau-t is the time between randomization - of all molecules. Inhibits particle dynamics somewhat, but little or + 'andersen-massive': Temperature coupling by randomizing velocities of all + particles at infrequent timesteps. Reference temperature and coupling + groups are selected as above. tau-t is the time between randomization + of all molecules. Inhibits particle dynamics somewhat, but little or no ergodicity issues. Currently only implemented with velocity Verlet. - 'v-rescale': Temperature coupling using velocity rescaling with a - stochastic term (JCP 126, 014101). This thermostat is similar to - Berendsen coupling, with the same scaling using tau-t, but the - stochastic term ensures that a proper canonical ensemble is generated. - The random seed is set with ld-seed. This thermostat works correctly - even for tau-t =0. For NVT simulations the conserved energy quantity + 'v-rescale': Temperature coupling using velocity rescaling with a + stochastic term (JCP 126, 014101). This thermostat is similar to + Berendsen coupling, with the same scaling using tau-t, but the + stochastic term ensures that a proper canonical ensemble is generated. + The random seed is set with ld-seed. This thermostat works correctly + even for tau-t =0. For NVT simulations the conserved energy quantity is written to the energy and log file. Default 'no' (since for the default integrator, 'sd', this option is ignored). """ nsttcouple: int = -1 """ - The frequency for coupling the temperature. The default value of -1 sets - nsttcouple equal to 100, or fewer steps if required for accurate + The frequency for coupling the temperature. The default value of -1 sets + nsttcouple equal to 100, or fewer steps if required for accurate integration (5 steps per tau for first order coupling, 20 steps per tau for - second order coupling). Note that the default value is large in order to - reduce the overhead of the additional computation and communication - required for obtaining the kinetic energy. For velocity Verlet integrators + second order coupling). Note that the default value is large in order to + reduce the overhead of the additional computation and communication + required for obtaining the kinetic energy. For velocity Verlet integrators nsttcouple is set to 1. Default -1. """ @@ -386,7 +385,7 @@ class Config: """ tau-t: FloatQuantity['picosecond'] = 2.0 * unit.picosecond """" - Time constant for coupling (one for each group in tc-grps), -1 means no + Time constant for coupling (one for each group in tc-grps), -1 means no temperature coupling. Default 2.0 * unit.picosecond. """ ref-t: FloatQuantity['Kelvin'] = 298.15 * unit.kelvin @@ -400,67 +399,67 @@ class Config: """ Options for pressure coupling (barostat). Allowed values are: 'no': No pressure coupling. This means a fixed box size. - 'berendsen': Exponential relaxation pressure coupling with time constant - tau-p. The box is scaled every nstpcouple steps. This barostat does not - yield a correct thermodynamic ensemble; it is only included to be able - to reproduce previous runs, and we strongly recommend against using it + 'berendsen': Exponential relaxation pressure coupling with time constant + tau-p. The box is scaled every nstpcouple steps. This barostat does not + yield a correct thermodynamic ensemble; it is only included to be able + to reproduce previous runs, and we strongly recommend against using it for new simulations. - 'C-rescale': Exponential relaxation pressure coupling with time constant - tau-p, including a stochastic term to enforce correct volume - fluctuations. The box is scaled every nstpcouple steps. It can be used + 'C-rescale': Exponential relaxation pressure coupling with time constant + tau-p, including a stochastic term to enforce correct volume + fluctuations. The box is scaled every nstpcouple steps. It can be used for both equilibration and production. - 'Parrinello-Rahman': Extended-ensemble pressure coupling where the box - vectors are subject to an equation of motion. The equation of motion - for the atoms is coupled to this. No instantaneous scaling takes place. - As for Nose-Hoover temperature coupling the time constant tau-p is the + 'Parrinello-Rahman': Extended-ensemble pressure coupling where the box + vectors are subject to an equation of motion. The equation of motion + for the atoms is coupled to this. No instantaneous scaling takes place. + As for Nose-Hoover temperature coupling the time constant tau-p is the period of pressure fluctuations at equilibrium. - 'MTTK': Martyna-Tuckerman-Tobias-Klein implementation, only useable with - integrator=md-vv or integrator=md-vv-avek, very similar to Parrinello-Rahman. - As for Nose-Hoover temperature coupling the time constant tau-p is the - period of pressure fluctuations at equilibrium. This is probably a - better method when you want to apply pressure scaling during data - collection, but beware that you can get very large oscillations if you - are starting from a different pressure. This requires a constant ensemble - temperature for the system. Currently it only supports isotropic scaling, + 'MTTK': Martyna-Tuckerman-Tobias-Klein implementation, only useable with + integrator=md-vv or integrator=md-vv-avek, very similar to Parrinello-Rahman. + As for Nose-Hoover temperature coupling the time constant tau-p is the + period of pressure fluctuations at equilibrium. This is probably a + better method when you want to apply pressure scaling during data + collection, but beware that you can get very large oscillations if you + are starting from a different pressure. This requires a constant ensemble + temperature for the system. Currently it only supports isotropic scaling, and only works without constraints. Default 'no'. """ pcoupltype: Literal['isotropic', 'semiisotropic', 'anisotropic', 'surface-tension'] = 'isotropic' """ - Specifies the kind of isotropy of the pressure coupling used. Each kind - takes one or more values for compressibility and ref-p. Only a single value + Specifies the kind of isotropy of the pressure coupling used. Each kind + takes one or more values for compressibility and ref-p. Only a single value is permitted for tau-p. Allowed values are: - 'isotropic': sotropic pressure coupling with time constant tau-p. One value + 'isotropic': sotropic pressure coupling with time constant tau-p. One value each for compressibility and ref-p is required. - 'semiisotropic': Pressure coupling which is isotropic in the x and y - direction, but different in the z direction. This can be useful for - membrane simulations. Two values each for compressibility and ref-p are + 'semiisotropic': Pressure coupling which is isotropic in the x and y + direction, but different in the z direction. This can be useful for + membrane simulations. Two values each for compressibility and ref-p are required, for x/y and z directions respectively. - 'anisotropic': Same as before, but 6 values are needed for xx, yy, zz, - xy/yx, xz/zx and yz/zy components, respectively. When the off-diagonal - compressibilities are set to zero, a rectangular box will stay - rectangular. Beware that anisotropic scaling can lead to extreme + 'anisotropic': Same as before, but 6 values are needed for xx, yy, zz, + xy/yx, xz/zx and yz/zy components, respectively. When the off-diagonal + compressibilities are set to zero, a rectangular box will stay + rectangular. Beware that anisotropic scaling can lead to extreme deformation of the simulation box. - 'surface-tension': Surface tension coupling for surfaces parallel to the - xy-plane. Uses normal pressure coupling for the z-direction, while the - surface tension is coupled to the x/y dimensions of the box. The first - ref-p value is the reference surface tension times the number of - surfaces bar nm, the second value is the reference z-pressure bar. - The two compressibility values are the compressibility in the x/y and z - direction respectively. The value for the z-compressibility should be - reasonably accurate since it influences the convergence of the - surface-tension, it can also be set to zero to have a box with constant + 'surface-tension': Surface tension coupling for surfaces parallel to the + xy-plane. Uses normal pressure coupling for the z-direction, while the + surface tension is coupled to the x/y dimensions of the box. The first + ref-p value is the reference surface tension times the number of + surfaces bar nm, the second value is the reference z-pressure bar. + The two compressibility values are the compressibility in the x/y and z + direction respectively. The value for the z-compressibility should be + reasonably accurate since it influences the convergence of the + surface-tension, it can also be set to zero to have a box with constant height. Default 'isotropic' """ nstpcouple: int = -1 """ - The frequency for coupling the pressure. The default value of -1 sets - nstpcouple equal to 100, or fewer steps if required for accurate integration - (5 steps per tau for first order coupling, 20 steps per tau for second order - coupling). Note that the default value is large in order to reduce the - overhead of the additional computation and communication required for - obtaining the virial and kinetic energy. For velocity Verlet integrators + The frequency for coupling the pressure. The default value of -1 sets + nstpcouple equal to 100, or fewer steps if required for accurate integration + (5 steps per tau for first order coupling, 20 steps per tau for second order + coupling). Note that the default value is large in order to reduce the + overhead of the additional computation and communication required for + obtaining the virial and kinetic energy. For velocity Verlet integrators nsttcouple is set to 1. Default -1. """ tau-p: FloatQuantity['picosecond'] = 5 * unit.picosecond @@ -470,26 +469,26 @@ class Config: """ compressibility: FloatQuantity['1/bar'] = 4.5e-05 / unit.bar """ - The compressibility. For water at 1 atm and 300 K the compressibility is + The compressibility. For water at 1 atm and 300 K the compressibility is 4.5e-5 bar-1. The number of required values is implied by pcoupltype. Default 4.5e-05 / unit.bar """ ref-p: FloatQuantity['bar'] = 1.01325 * unit.bar """ - The reference pressure for coupling. The number of required values is + The reference pressure for coupling. The number of required values is implied by pcoupltype. Default 1.01325 * unit.bar. """ refcoord-scaling: Literal['no', 'all', 'com'] = 'no' """ Allowed values are: 'no': The reference coordinates for position restraints are not modified. - 'all': The reference coordinates are scaled with the scaling matrix of the + 'all': The reference coordinates are scaled with the scaling matrix of the pressure coupling. - 'com': Scale the center of mass of the reference coordinates with the - scaling matrix of the pressure coupling. The vectors of each reference - coordinate to the center of mass are not scaled. Only one COM is used, - even when there are multiple molecules with position restraints. - For calculating the COM of the reference coordinates in the starting + 'com': Scale the center of mass of the reference coordinates with the + scaling matrix of the pressure coupling. The vectors of each reference + coordinate to the center of mass are not scaled. Only one COM is used, + even when there are multiple molecules with position restraints. + For calculating the COM of the reference coordinates in the starting configuration, periodic boundary conditions are not taken into account. Default 'no'. """ @@ -498,10 +497,10 @@ class Config: gen-vel: Literal['no', 'yes'] = 'yes' """ Velocity generation. Allowed values are: - 'no': Do not generate velocities. The velocities are set to zero when there + 'no': Do not generate velocities. The velocities are set to zero when there are no velocities in the input structure file. - 'yes': Generate velocities in gmx grompp according to a Maxwell distribution - at temperature gen-temp, with random seed gen-seed. This is only + 'yes': Generate velocities in gmx grompp according to a Maxwell distribution + at temperature gen-temp, with random seed gen-seed. This is only meaningful with integrator=md. Default 'yes'. """ @@ -511,21 +510,21 @@ class Config: """ gen-seed: int = -1 """ - Used to initialize random generator for random velocities, when gen-seed is + Used to initialize random generator for random velocities, when gen-seed is set to -1, a pseudo random seed is used. Default -1. """ ### Bonds ### constraints: Literal['none', 'h-bonds', 'all-bonds', 'h-angles', 'all-angles'] = 'h-bonds' """ - Controls which bonds in the topology will be converted to rigid holonomic - constraints. Note that typical rigid water models do not have bonds, but + Controls which bonds in the topology will be converted to rigid holonomic + constraints. Note that typical rigid water models do not have bonds, but rather a specialized [settles] directive, so are not affected by this keyword. Allowed values are: 'none': No bonds converted to constraints. 'h-bonds': Convert the bonds with H-atoms to constraints. 'all-bonds': Convert all bonds to constraints. - 'h-angles': Convert all bonds to constraints and convert the angles that + 'h-angles': Convert all bonds to constraints and convert the angles that involve H-atoms to bond-constraints. 'all-angles': Convert all bonds to constraints and all angles to bond-constraints. @@ -535,22 +534,22 @@ class Config: """ Chooses which solver satisfies any non-SETTLE holonomic constraints. Allowed values are: - 'lincs': LINear Constraint Solver. With domain decomposition the parallel - version P-LINCS is used. The accuracy in set with lincs-order, which - sets the number of matrices in the expansion for the matrix inversion. - After the matrix inversion correction the algorithm does an iterative - correction to compensate for lengthening due to rotation. The number of - such iterations can be controlled with lincs-iter. The root mean square - relative constraint deviation is printed to the log file every nstlog - steps. If a bond rotates more than lincs-warnangle in one step, a - warning will be printed both to the log file and to stderr. + 'lincs': LINear Constraint Solver. With domain decomposition the parallel + version P-LINCS is used. The accuracy in set with lincs-order, which + sets the number of matrices in the expansion for the matrix inversion. + After the matrix inversion correction the algorithm does an iterative + correction to compensate for lengthening due to rotation. The number of + such iterations can be controlled with lincs-iter. The root mean square + relative constraint deviation is printed to the log file every nstlog + steps. If a bond rotates more than lincs-warnangle in one step, a + warning will be printed both to the log file and to stderr. LINCS should not be used with coupled angle constraints. - 'shake': SHAKE is slightly slower and less stable than LINCS, but does work - with angle constraints. The relative tolerance is set with shake-tol, - 0.0001 is a good value for “normal” MD. SHAKE does not support - constraints between atoms on different decomposition domains, so it can - only be used with domain decomposition when so-called update-groups are - used, which is usally the case when only bonds involving hydrogens are + 'shake': SHAKE is slightly slower and less stable than LINCS, but does work + with angle constraints. The relative tolerance is set with shake-tol, + 0.0001 is a good value for “normal” MD. SHAKE does not support + constraints between atoms on different decomposition domains, so it can + only be used with domain decomposition when so-called update-groups are + used, which is usally the case when only bonds involving hydrogens are constrained. SHAKE can not be used with energy minimization. Default 'lincs' """ @@ -559,25 +558,25 @@ class Config: This option was formerly known as unconstrained-start. Allowed values are: 'no': Apply constraints to the start configuration and reset shells. - 'yes': Do not apply constraints to the start configuration and do not reset + 'yes': Do not apply constraints to the start configuration and do not reset shells, useful for exact coninuation and reruns. """ shake-tol: float = 0.0001 """ Relative tolerance for SHAKE. Default 0.0001. """ lincs-order: int = 12 """ - Highest order in the expansion of the constraint coupling matrix. When - constraints form triangles, an additional expansion of the same order is - applied on top of the normal expansion only for the couplings within such - triangles. For “normal” MD simulations an order of 4 usually suffices, 6 is - needed for large time-steps with virtual sites or BD. For accurate energy - minimization in double precision an order of 8 or more might be required. - Note that in single precision an order higher than 6 will often lead to + Highest order in the expansion of the constraint coupling matrix. When + constraints form triangles, an additional expansion of the same order is + applied on top of the normal expansion only for the couplings within such + triangles. For “normal” MD simulations an order of 4 usually suffices, 6 is + needed for large time-steps with virtual sites or BD. For accurate energy + minimization in double precision an order of 8 or more might be required. + Note that in single precision an order higher than 6 will often lead to worse accuracy due to amplification of rounding errors. Default 12. """ lincs-iter: int = 1 """ - Number of iterations to correct for rotational lengthening in LINCS. + Number of iterations to correct for rotational lengthening in LINCS. Default 1. """ lincs-warnangle: FloatQuantity['deg'] = 30 * unit.degree @@ -589,7 +588,7 @@ class Config: """ Allowed options are: 'no': Bonds are represented by a harmonic potential. - 'yes': Bonds are represented by a Morse potential. + 'yes': Bonds are represented by a Morse potential. Default 'no'. """ @@ -599,50 +598,50 @@ class MDOutputSettings(SettingsBaseModel): """ nstxout: int = 0 """ - Number of steps that elapse between writing coordinates to the output - trajectory file (trr), the last coordinates are always written unless 0, + Number of steps that elapse between writing coordinates to the output + trajectory file (trr), the last coordinates are always written unless 0, which means coordinates are not written into the trajectory file. Default 0. """ nstvout: int = 0 """ - Number of steps that elapse between writing velocities to the output - trajectory file (trr), the last velocities are always written unless 0, + Number of steps that elapse between writing velocities to the output + trajectory file (trr), the last velocities are always written unless 0, which means velocities are not written into the trajectory file. Default 0. """ nstfout: int = 0 """ Number of steps that elapse between writing forces to the output trajectory - file (trr), the last forces are always written, unless 0, which means + file (trr), the last forces are always written, unless 0, which means forces are not written into the trajectory file. Default 0. """ nstlog: int = 1000 """ - Number of steps that elapse between writing energies to the log file, the + Number of steps that elapse between writing energies to the log file, the last energies are always written. Default 1000. """ nstcalcenergy: int = 100 """ - Number of steps that elapse between calculating the energies, 0 is never. - This option is only relevant with dynamics. This option affects the - performance in parallel simulations, because calculating energies requires + Number of steps that elapse between calculating the energies, 0 is never. + This option is only relevant with dynamics. This option affects the + performance in parallel simulations, because calculating energies requires global communication between all processes which can become a bottleneck at high parallelization. Default 100. """ nstenergy: int = 1000 """" - Number of steps that elapse between writing energies to energy file, the - last energies are always written, should be a multiple of nstcalcenergy. - Note that the exact sums and fluctuations over all MD steps modulo - nstcalcenergy are stored in the energy file, so gmx energy can report + Number of steps that elapse between writing energies to energy file, the + last energies are always written, should be a multiple of nstcalcenergy. + Note that the exact sums and fluctuations over all MD steps modulo + nstcalcenergy are stored in the energy file, so gmx energy can report exact energy averages and fluctuations also when nstenergy > 1 """ nstxout-compressed: int = 0 """ - Number of steps that elapse between writing position coordinates using - lossy compression (xtc file), 0 for not writing compressed coordinates + Number of steps that elapse between writing position coordinates using + lossy compression (xtc file), 0 for not writing compressed coordinates output. Default 0. """ compressed-x-precision: int = 1000 @@ -652,12 +651,12 @@ class MDOutputSettings(SettingsBaseModel): """ compressed-x-grps: str """ - Group(s) to write to the compressed trajectory file, by default the whole + Group(s) to write to the compressed trajectory file, by default the whole system is written (if nstxout-compressed > 0). """ energygrps: str """ - Group(s) for which to write to write short-ranged non-bonded potential + Group(s) for which to write to write short-ranged non-bonded potential energies to the energy file (not supported on GPUs) """ From 59fe6d92aa93cda6c81c7fe323fa9f5cc0bea4cd Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Mon, 15 Jul 2024 15:21:31 +0200 Subject: [PATCH 04/23] Separate settings for EM, NVT, NPT --- .../protocols/gromacs_md/md_settings.py | 86 +++++++++++++------ 1 file changed, 61 insertions(+), 25 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 708cb54..bf8aa6e 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -26,14 +26,14 @@ except ImportError: from pydantic import validator # type: ignore[assignment] -class MDSimulationSettings(SettingsBaseModel): +class SimulationSettings(SettingsBaseModel): """ - Settings for MD simulations in Gromacs. + Settings for simulations in Gromacs. """ class Config: arbitrary_types_allowed = True - ### Run control ### + ###Run control### integrator: Literal['md', 'md-vv', 'md-vv-avek', 'sd', 'steep'] = 'sd' """ MD integrators and other algorithms (steep for energy minimization). @@ -105,7 +105,7 @@ class Config: Group(s) for center of mass motion removal, default is the whole system. """ - ### Langevin dynamics ### + ###Langevin dynamics### ld-seed: int = -1 """ Integer used to initialize random generator for thermal noise for @@ -113,17 +113,6 @@ class Config: a pseudo random seed is used. Default -1. """ - ### Energy minimization ### - emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) - """ - The minimization is converged when the maximum force is smaller than this - value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) - """ - emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer - """ - Initial step size. Default 0.01 * unit.nanometer - """ - ### Neighbor searching ### cutoff-scheme: Literal['verlet'] = 'verlet' """ @@ -173,7 +162,7 @@ class Config: options and the value of rlist is ignored. """ - ### Electrostatics ## + ###Electrostatics### coulombtype: Literal['cut-off', 'ewald', 'PME', 'P3M-AD', 'Reaction-Field'] = 'PME' """ Treatment of electrostatics @@ -228,7 +217,7 @@ class Config: with reaction-field electrostatics. A value of 0 means infinity. Default 0 """ - ### Van der Waals ### + ###Van der Waals### vdwtype: Literal['Cut-off', 'PME'] = 'Cut-off' """ Treatment of vdW interactions. Allowed options are: @@ -274,7 +263,7 @@ class Config: Default 'EnerPres' """ - ### Ewald ### + ###Ewald### fourierspacing: FloatQuantity['nanometer'] = 0.12 * unit.nanometer """ For ordinary Ewald, the ratio of the box dimensions and the spacing @@ -330,7 +319,7 @@ class Config: variant of the long range corrections. Default 0. """ - ### Temerature coupling ### + ###Temperature coupling### tcoupl: Literal['no', 'berendsen', 'nose-hoover', 'andersen', 'andersen-massive', 'v-rescale'] = 'no' """ Temperature coupling options. Note that tcoupl will be ignored when the @@ -395,7 +384,7 @@ class Config: Default 298.15 * unit.kelvin """ - ### Pressure coupling ### + ###Pressure coupling### pcoupl: Literal['no', 'berendsen', 'C-rescale', 'Parrinello-Rahman', 'MTTK'] = 'no' """ Options for pressure coupling (barostat). Allowed values are: @@ -494,7 +483,7 @@ class Config: Default 'no'. """ - ### Velocity generation ### + ###Velocity generation### gen-vel: Literal['no', 'yes'] = 'yes' """ Velocity generation. Allowed values are: @@ -515,7 +504,7 @@ class Config: set to -1, a pseudo random seed is used. Default -1. """ - ### Bonds ### + ###Bonds### constraints: Literal['none', 'h-bonds', 'all-bonds', 'h-angles', 'all-angles'] = 'h-bonds' """ Controls which bonds in the topology will be converted to rigid holonomic @@ -593,7 +582,7 @@ class Config: Default 'no'. """ -class MDOutputSettings(SettingsBaseModel): +class OutputSettings(SettingsBaseModel): """" Output Settings for simulations run using Gromacs """ @@ -661,6 +650,49 @@ class MDOutputSettings(SettingsBaseModel): energies to the energy file (not supported on GPUs) """ +class EMSimulationSettings(SimulationSettings): + """ + Settings for energy minimization. + """ + integrator = 'steep' + nsteps = 5000 + + emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) + """ + The minimization is converged when the maximum force is smaller than this + value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) + """ + emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer + """ + Initial step size. Default 0.01 * unit.nanometer + """ + +class EMOutputSettings(OutputSettings): + """ + Output Settings for the energy minimization. + """ + +class NVTSimulationSettings(SimulationSettings): + """ + Settings for MD simulation in the NVT ensemble. + """ + +class NVTOutputSettings(OutputSettings): + """ + Output Settings for the MD simulation in the NVT ensemble. + """ + +class NPTSimulationSettings(SimulationSettings): + """ + Settings for MD simulation in the NPT ensemble. + """ + +class NPTOutputSettings(OutputSettings): + """ + Output Settings for the MD simulation in the NPT ensemble. + """ + + class GromacsMDProtocolSettings(Settings): class Config: @@ -690,7 +722,11 @@ def must_be_positive(cls, v): integrator_settings: IntegratorSettings # Simulation run settings - simulation_settings: MDSimulationSettings + simulation_settings_em: EMSimulationSettings + simulation_settings_nvt: NVTSimulationSettings + simulation_settings_npt: NPTSimulationSettings # Simulations output settings - output_settings: MDOutputSettings + output_settings_em: EMOutputSettings + output_settings_nvt: NVTOutputSettings + output_settings_npt: NPTOutputSettings From 185628aea7b14bd4832824aced3751c4c480fbaf Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 15 Jul 2024 13:37:21 +0000 Subject: [PATCH 05/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 5d77682..cd639c3 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -672,7 +672,7 @@ class EMSimulationSettings(SimulationSettings): emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) """ - The minimization is converged when the maximum force is smaller than this + The minimization is converged when the maximum force is smaller than this value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) """ emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer From 797fe5df5b9596778e06aba8e350781a59700013 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Mon, 15 Jul 2024 15:53:51 +0200 Subject: [PATCH 06/23] fix pep8 issues --- .../protocols/gromacs_md/md_settings.py | 173 +++++++++--------- 1 file changed, 84 insertions(+), 89 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index cd639c3..0f134b3 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -25,6 +25,7 @@ except ImportError: from pydantic import validator # type: ignore[assignment] + class SimulationSettings(SettingsBaseModel): """ Settings for simulations in Gromacs. @@ -32,7 +33,7 @@ class SimulationSettings(SettingsBaseModel): class Config: arbitrary_types_allowed = True - ###Run control### + # # # Run control # # # integrator: Literal['md', 'md-vv', 'md-vv-avek', 'sd', 'steep'] = 'sd' """ MD integrators and other algorithms (steep for energy minimization). @@ -63,13 +64,13 @@ class Config: Maximum number of steps to integrate or minimize, -1 is no maximum Default 0 """ - init-step: int = 0 + init_step: int = 0 """ The starting step. The time at step i in a run is calculated as: t = tinit + dt * (init-step + i). Default 0 """ - simulation-part: int=0 + simulation_part: int = 0 """ A simulation can consist of multiple parts, each of which has a part number. This option specifies what that number will be, which helps keep track of @@ -77,13 +78,13 @@ class Config: useful to set only when coping with a crashed simulation where files were lost. Default 0 """ - mass-repartition-factor: int = 1 + mass_repartition_factor: int = 1 """ Scales the masses of the lightest atoms in the system by this factor to the mass mMin. All atoms with a mass lower than mMin also have their mass set to that mMin. Default 1 (no mass scaling) """ - comm-mode: Literal['Linear', 'Angular', 'Linear-acceleration-correction', 'None'] = 'Linear' + comm_mode: Literal['Linear', 'Angular', 'Linear-acceleration-correction', 'None'] = 'Linear' """ Settings for center of mass treatmeant Allowed values are: @@ -99,35 +100,21 @@ class Config: Frequency for center of mass motion removal in unit [steps]. Default 100 """ - comm-grps: str + comm_grps: str """ Group(s) for center of mass motion removal, default is the whole system. """ - ###Langevin dynamics### - ld-seed: int = -1 + # # # Langevin dynamics # # # + ld_seed: int = -1 """ Integer used to initialize random generator for thermal noise for stochastic and Brownian dynamics. When ld-seed is set to -1, a pseudo random seed is used. Default -1. """ -<<<<<<< HEAD -======= - ### Energy minimization ### - emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) - """ - The minimization is converged when the maximum force is smaller than this - value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) - """ - emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer - """ - Initial step size. Default 0.01 * unit.nanometer - """ - ->>>>>>> 10d6f0963fe5c092695dbee628812d518a30d37a - ### Neighbor searching ### - cutoff-scheme: Literal['verlet'] = 'verlet' + # # # Neighbor searching # # # + cutoff_scheme: Literal['verlet'] = 'verlet' """ Only allowed option: 'verlet': Generate a pair list with buffering. The buffer size is @@ -152,14 +139,14 @@ class Config: This can be used in combination with walls. Default 'xyz'. """ - verlet-buffer-tolerance: FloatQuantity['kilojoule / (mole * picosecond)'] = 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) + verlet_buffer_tolerance: FloatQuantity['kilojoule / (mole * picosecond)'] = 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) """" Used when performing a simulation with dynamics. This sets the maximum allowed error for pair interactions per particle caused by the Verlet buffer, which indirectly sets rlist. Default 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) """ - verlet-buffer-pressure-tolerance: FloatQuantity['bar'] = 0.5 * unit.bar + verlet_buffer_pressure_tolerance: FloatQuantity['bar'] = 0.5 * unit.bar """ Used when performing a simulation with dynamics and only active when verlet-buffer-tolerance is positive. This sets the maximum tolerated error @@ -171,11 +158,11 @@ class Config: rlist: FloatQuantity['nanometer'] = 1 * unit.nanometer """ Cut-off distance for the short-range neighbor list. With dynamics, this is - by default set by the verlet-buffer-tolerance and verlet-buffer-pressure-tolerance - options and the value of rlist is ignored. + by default set by the verlet-buffer-tolerance and + verlet-buffer-pressure-tolerance options and the value of rlist is ignored. """ - ###Electrostatics### + # # # Electrostatics # # # coulombtype: Literal['cut-off', 'ewald', 'PME', 'P3M-AD', 'Reaction-Field'] = 'PME' """ Treatment of electrostatics @@ -199,7 +186,7 @@ class Config: infinity by setting epsilon-rf =0. Default 'PME' """ - coulomb-modifier: Literal['Potential-shift', 'None'] = 'Potential-shift' + coulomb_modifier: Literal['Potential-shift', 'None'] = 'Potential-shift' """ Allowed options are: 'Potential-shift': Shift the Coulomb potential by a constant such that it @@ -208,7 +195,7 @@ class Config: comparing energies with those computed with other software. Default 'Potential-shift' """ - rcoulomb-switch: FloatQuantity['nanometer'] = 0 * unit.nanometer + rcoulomb_switch: FloatQuantity['nanometer'] = 0 * unit.nanometer """ Where to start switching the Coulomb potential, only relevant when force or potential switching is used. @@ -220,17 +207,17 @@ class Config: increased by the PME tuning in gmx mdrun along with the PME grid spacing. Default 1.2 * unit.nanometer """ - epsilon-r: int = 1 + epsilon_r: int = 1 """ The relative dielectric constant. A value of 0 means infinity. Default 1. """ - epsilon-rf: int = 0 + epsilon_rf: int = 0 """ The relative dielectric constant of the reaction field. This is only used with reaction-field electrostatics. A value of 0 means infinity. Default 0 """ - ###Van der Waals### + # # # Van der Waals # # # vdwtype: Literal['Cut-off', 'PME'] = 'Cut-off' """ Treatment of vdW interactions. Allowed options are: @@ -239,7 +226,7 @@ class Config: 'PME: Fast smooth Particle-mesh Ewald (SPME) for VdW interactions. Default 'Cut-off'. """ - vdw-modifier: Literal['Potential-shift', 'None', 'Force-switch', 'Potential-switch'] = 'Potential-shift' + vdw_modifier: Literal['Potential-shift', 'None', 'Force-switch', 'Potential-switch'] = 'Potential-shift' """ Allowed values are: 'Potential-shift': Shift the Van der Waals potential by a constant such @@ -257,7 +244,7 @@ class Config: This option should only be used if the force field you are using requires this. """ - rvdw-switch: FloatQuantity['nanometer'] = 0 * unit.nanometer + rvdw_switch: FloatQuantity['nanometer'] = 0 * unit.nanometer """" Where to start switching the LJ force and possibly the potential, only relevant when force or potential switching is used. @@ -276,7 +263,7 @@ class Config: Default 'EnerPres' """ - ###Ewald### + # # # Ewald # # # fourierspacing: FloatQuantity['nanometer'] = 0.12 * unit.nanometer """ For ordinary Ewald, the ratio of the box dimensions and the spacing @@ -287,24 +274,24 @@ class Config: entering a non-zero value for that fourier-nx direction. Default 0.12 * unit.nanometer """ - pme-order: int = 4 + pme_order: int = 4 """ The number of grid points along a dimension to which a charge is mapped. The actual order of the PME interpolation is one less, e.g. the default of - 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for P3M-AD). - When running in parallel, it can be worth to switch to 5 and simultaneously - increase the grid spacing. Note that on the CPU only values 4 and 5 have - SIMD acceleration and GPUs only support the value 4. + 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for + P3M-AD). When running in parallel, it can be worth to switch to 5 and + simultaneously increase the grid spacing. Note that on the CPU only values + 4 and 5 have SIMD acceleration and GPUs only support the value 4. Default 4. """ - ewald-rtol: float = 1e-5 + ewald_rtol: float = 1e-5 """ The relative strength of the Ewald-shifted direct potential at rcoulomb is given by ewald-rtol. Decreasing this will give a more accurate direct sum, but then you need more wave vectors for the reciprocal sum. Default 1e-5 """ - lj-pme-comb-rule: Literal['Geometric', 'Lorentz-Berthelot'] = 'Geometric' + lj_pme_comb_rule: Literal['Geometric', 'Lorentz-Berthelot'] = 'Geometric' """ The combination rules used to combine VdW-parameters in the reciprocal part of LJ-PME. Geometric rules are much faster than Lorentz-Berthelot and @@ -314,15 +301,16 @@ class Config: 'Lorentz-Berthelot': Apply Lorentz-Berthelot combination rules. Default 'Geometric'. """ - ewald-geometry: Literal['3d', '3dc'] = '3d' + ewald_geometry: Literal['3d', '3dc'] = '3d' """ Allowed values are: '3d': The Ewald sum is performed in all three dimensions. '3dc': The reciprocal sum is still performed in 3D, but a force and - potential correction applied in the z dimension to produce a pseudo-2D summation. + potential correction applied in the z dimension to produce a pseudo-2D + summation. Default '3d'. """ - epsilon-surface: float = 0 + epsilon_surface: float = 0 """ This controls the dipole correction to the Ewald summation in 3D. The default value of zero means it is turned off. Turn it on by setting it @@ -332,7 +320,7 @@ class Config: variant of the long range corrections. Default 0. """ - ###Temperature coupling### + # # # Temperature coupling # # # tcoupl: Literal['no', 'berendsen', 'nose-hoover', 'andersen', 'andersen-massive', 'v-rescale'] = 'no' """ Temperature coupling options. Note that tcoupl will be ignored when the @@ -369,7 +357,7 @@ class Config: The random seed is set with ld-seed. This thermostat works correctly even for tau-t =0. For NVT simulations the conserved energy quantity is written to the energy and log file. - Default 'no' (since for the default integrator, 'sd', this option is ignored). + Default 'no' (for the default integrator, 'sd', this option is ignored). """ nsttcouple: int = -1 """ @@ -382,22 +370,22 @@ class Config: nsttcouple is set to 1. Default -1. """ - tc-grps: str = 'system' + tc_grps: str = 'system' """ Groups to couple to separate temperature baths. Default 'system'. """ - tau-t: FloatQuantity['picosecond'] = 2.0 * unit.picosecond + tau_t: FloatQuantity['picosecond'] = 2.0 * unit.picosecond """" Time constant for coupling (one for each group in tc-grps), -1 means no temperature coupling. Default 2.0 * unit.picosecond. """ - ref-t: FloatQuantity['Kelvin'] = 298.15 * unit.kelvin + ref_t: FloatQuantity['Kelvin'] = 298.15 * unit.kelvin """ Reference temperature for coupling (one for each group in tc-grps). Default 298.15 * unit.kelvin """ - ###Pressure coupling### + # # # Pressure coupling # # # pcoupl: Literal['no', 'berendsen', 'C-rescale', 'Parrinello-Rahman', 'MTTK'] = 'no' """ Options for pressure coupling (barostat). Allowed values are: @@ -417,14 +405,15 @@ class Config: As for Nose-Hoover temperature coupling the time constant tau-p is the period of pressure fluctuations at equilibrium. 'MTTK': Martyna-Tuckerman-Tobias-Klein implementation, only useable with - integrator=md-vv or integrator=md-vv-avek, very similar to Parrinello-Rahman. - As for Nose-Hoover temperature coupling the time constant tau-p is the - period of pressure fluctuations at equilibrium. This is probably a - better method when you want to apply pressure scaling during data - collection, but beware that you can get very large oscillations if you - are starting from a different pressure. This requires a constant ensemble - temperature for the system. Currently it only supports isotropic scaling, - and only works without constraints. + integrator=md-vv or integrator=md-vv-avek, very similar to + Parrinello-Rahman. As for Nose-Hoover temperature coupling the time + constant tau-p is the period of pressure fluctuations at equilibrium. + This is probably a better method when you want to apply pressure + scaling during data collection, but beware that you can get very large + oscillations if you are starting from a different pressure. + This requires a constant ensemble temperature for the system. + Currently it only supports isotropic scaling, and only works without + constraints. Default 'no'. """ pcoupltype: Literal['isotropic', 'semiisotropic', 'anisotropic', 'surface-tension'] = 'isotropic' @@ -458,14 +447,14 @@ class Config: nstpcouple: int = -1 """ The frequency for coupling the pressure. The default value of -1 sets - nstpcouple equal to 100, or fewer steps if required for accurate integration - (5 steps per tau for first order coupling, 20 steps per tau for second order - coupling). Note that the default value is large in order to reduce the - overhead of the additional computation and communication required for - obtaining the virial and kinetic energy. For velocity Verlet integrators - nsttcouple is set to 1. Default -1. + nstpcouple equal to 100, or fewer steps if required for accurate + integration (5 steps per tau for first order coupling, 20 steps per tau for + second order coupling). Note that the default value is large in order to + reduce the overhead of the additional computation and communication + required for obtaining the virial and kinetic energy. For velocity Verlet + integrators nsttcouple is set to 1. Default -1. """ - tau-p: FloatQuantity['picosecond'] = 5 * unit.picosecond + tau_p: FloatQuantity['picosecond'] = 5 * unit.picosecond """ The time constant for pressure coupling (one value for all directions). Default 5 * unit.picosecond. @@ -476,12 +465,12 @@ class Config: 4.5e-5 bar-1. The number of required values is implied by pcoupltype. Default 4.5e-05 / unit.bar """ - ref-p: FloatQuantity['bar'] = 1.01325 * unit.bar + ref_p: FloatQuantity['bar'] = 1.01325 * unit.bar """ The reference pressure for coupling. The number of required values is implied by pcoupltype. Default 1.01325 * unit.bar. """ - refcoord-scaling: Literal['no', 'all', 'com'] = 'no' + refcoord_scaling: Literal['no', 'all', 'com'] = 'no' """ Allowed values are: 'no': The reference coordinates for position restraints are not modified. @@ -496,34 +485,34 @@ class Config: Default 'no'. """ - ###Velocity generation### - gen-vel: Literal['no', 'yes'] = 'yes' + # # # Velocity generation # # # + gen_vel: Literal['no', 'yes'] = 'yes' """ Velocity generation. Allowed values are: 'no': Do not generate velocities. The velocities are set to zero when there are no velocities in the input structure file. - 'yes': Generate velocities in gmx grompp according to a Maxwell distribution - at temperature gen-temp, with random seed gen-seed. This is only - meaningful with integrator=md. + 'yes': Generate velocities in gmx grompp according to a Maxwell + distribution at temperature gen-temp, with random seed gen-seed. + This is only meaningful with integrator=md. Default 'yes'. """ - gen-temp: FloatQuantity['kelvin'] = 298.15 * unit.kelvin + gen_temp: FloatQuantity['kelvin'] = 298.15 * unit.kelvin """ Temperature for Maxwell distribution. Default 298.15 * unit.kelvin. """ - gen-seed: int = -1 + gen_seed: int = -1 """ Used to initialize random generator for random velocities, when gen-seed is set to -1, a pseudo random seed is used. Default -1. """ - ###Bonds### + # # # Bonds # # # constraints: Literal['none', 'h-bonds', 'all-bonds', 'h-angles', 'all-angles'] = 'h-bonds' """ Controls which bonds in the topology will be converted to rigid holonomic constraints. Note that typical rigid water models do not have bonds, but - rather a specialized [settles] directive, so are not affected by this keyword. - Allowed values are: + rather a specialized [settles] directive, so are not affected by this + keyword. Allowed values are: 'none': No bonds converted to constraints. 'h-bonds': Convert the bonds with H-atoms to constraints. 'all-bonds': Convert all bonds to constraints. @@ -533,7 +522,7 @@ class Config: bond-constraints. Default 'h-bonds' """ - constraint-algorithm: Literal['lincs', 'shake'] = 'lincs' + constraint_algorithm: Literal['lincs', 'shake'] = 'lincs' """ Chooses which solver satisfies any non-SETTLE holonomic constraints. Allowed values are: @@ -564,9 +553,9 @@ class Config: 'yes': Do not apply constraints to the start configuration and do not reset shells, useful for exact coninuation and reruns. """ - shake-tol: float = 0.0001 + shake_tol: float = 0.0001 """ Relative tolerance for SHAKE. Default 0.0001. """ - lincs-order: int = 12 + lincs_order: int = 12 """ Highest order in the expansion of the constraint coupling matrix. When constraints form triangles, an additional expansion of the same order is @@ -577,12 +566,12 @@ class Config: Note that in single precision an order higher than 6 will often lead to worse accuracy due to amplification of rounding errors. Default 12. """ - lincs-iter: int = 1 + lincs_iter: int = 1 """ Number of iterations to correct for rotational lengthening in LINCS. Default 1. """ - lincs-warnangle: FloatQuantity['deg'] = 30 * unit.degree + lincs_warnangle: FloatQuantity['deg'] = 30 * unit.degree """ Maximum angle that a bond can rotate before LINCS will complain. Default 30 * unit.degree. @@ -595,6 +584,7 @@ class Config: Default 'no'. """ + class OutputSettings(SettingsBaseModel): """" Output Settings for simulations run using Gromacs @@ -641,18 +631,18 @@ class OutputSettings(SettingsBaseModel): nstcalcenergy are stored in the energy file, so gmx energy can report exact energy averages and fluctuations also when nstenergy > 1 """ - nstxout-compressed: int = 0 + nstxout_compressed: int = 0 """ Number of steps that elapse between writing position coordinates using lossy compression (xtc file), 0 for not writing compressed coordinates output. Default 0. """ - compressed-x-precision: int = 1000 + compressed_x_precision: int = 1000 """ Precision with which to write to the compressed trajectory file. Default 1000. """ - compressed-x-grps: str + compressed_x_grps: str """ Group(s) to write to the compressed trajectory file, by default the whole system is written (if nstxout-compressed > 0). @@ -663,6 +653,7 @@ class OutputSettings(SettingsBaseModel): energies to the energy file (not supported on GPUs) """ + class EMSimulationSettings(SimulationSettings): """ Settings for energy minimization. @@ -680,33 +671,37 @@ class EMSimulationSettings(SimulationSettings): Initial step size. Default 0.01 * unit.nanometer """ + class EMOutputSettings(OutputSettings): """ Output Settings for the energy minimization. """ + class NVTSimulationSettings(SimulationSettings): """ Settings for MD simulation in the NVT ensemble. """ + class NVTOutputSettings(OutputSettings): """ Output Settings for the MD simulation in the NVT ensemble. """ + class NPTSimulationSettings(SimulationSettings): """ Settings for MD simulation in the NPT ensemble. """ + class NPTOutputSettings(OutputSettings): """ Output Settings for the MD simulation in the NPT ensemble. """ - class GromacsMDProtocolSettings(Settings): class Config: arbitrary_types_allowed = True From 1adfab6d8f7c2673738e05accf5e8c30f9b0a6bb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 15 Jul 2024 13:54:00 +0000 Subject: [PATCH 07/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../protocols/gromacs_md/md_settings.py | 147 ++++++++++-------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 0f134b3..65c8372 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -30,11 +30,12 @@ class SimulationSettings(SettingsBaseModel): """ Settings for simulations in Gromacs. """ + class Config: arbitrary_types_allowed = True # # # Run control # # # - integrator: Literal['md', 'md-vv', 'md-vv-avek', 'sd', 'steep'] = 'sd' + integrator: Literal["md", "md-vv", "md-vv-avek", "sd", "steep"] = "sd" """ MD integrators and other algorithms (steep for energy minimization). Allowed values are: @@ -48,13 +49,13 @@ class Config: integrator the parameters tcoupl and nsttcouple are ignored. 'steep': A steepest descent algorithm for energy minimization. """ - tinit: FloatQuantity['picosecond'] = 0 * unit.picosecond + tinit: FloatQuantity["picosecond"] = 0 * unit.picosecond """ Starting time for the MD run. This only makes sense for time-based integrators. Default 0 * unit.picosecond """ - dt: FloatQuantity['picosecond'] = 0.001 * unit.picosecond + dt: FloatQuantity["picosecond"] = 0.001 * unit.picosecond """ Time step for integration (only makes sense for time-based integrators). Default 0.001 * unit.picosecond @@ -84,7 +85,9 @@ class Config: the mass mMin. All atoms with a mass lower than mMin also have their mass set to that mMin. Default 1 (no mass scaling) """ - comm_mode: Literal['Linear', 'Angular', 'Linear-acceleration-correction', 'None'] = 'Linear' + comm_mode: Literal[ + "Linear", "Angular", "Linear-acceleration-correction", "None" + ] = "Linear" """ Settings for center of mass treatmeant Allowed values are: @@ -114,7 +117,7 @@ class Config: """ # # # Neighbor searching # # # - cutoff_scheme: Literal['verlet'] = 'verlet' + cutoff_scheme: Literal["verlet"] = "verlet" """ Only allowed option: 'verlet': Generate a pair list with buffering. The buffer size is @@ -129,7 +132,7 @@ class Config: If set to zero the neighbor list is only constructed once and never updated. Default 10. """ - pbc: Literal['xyz', 'no', 'xy'] = 'xyz' + pbc: Literal["xyz", "no", "xy"] = "xyz" """ Treatment of periodic boundary conditions. Allowed values are: @@ -139,14 +142,16 @@ class Config: This can be used in combination with walls. Default 'xyz'. """ - verlet_buffer_tolerance: FloatQuantity['kilojoule / (mole * picosecond)'] = 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) + verlet_buffer_tolerance: FloatQuantity["kilojoule / (mole * picosecond)"] = ( + 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) + ) """" Used when performing a simulation with dynamics. This sets the maximum allowed error for pair interactions per particle caused by the Verlet buffer, which indirectly sets rlist. Default 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) """ - verlet_buffer_pressure_tolerance: FloatQuantity['bar'] = 0.5 * unit.bar + verlet_buffer_pressure_tolerance: FloatQuantity["bar"] = 0.5 * unit.bar """ Used when performing a simulation with dynamics and only active when verlet-buffer-tolerance is positive. This sets the maximum tolerated error @@ -155,15 +160,15 @@ class Config: as the pair list ages. Default 0.5 * unit.bar """ - rlist: FloatQuantity['nanometer'] = 1 * unit.nanometer + rlist: FloatQuantity["nanometer"] = 1 * unit.nanometer """ Cut-off distance for the short-range neighbor list. With dynamics, this is - by default set by the verlet-buffer-tolerance and + by default set by the verlet-buffer-tolerance and verlet-buffer-pressure-tolerance options and the value of rlist is ignored. """ # # # Electrostatics # # # - coulombtype: Literal['cut-off', 'ewald', 'PME', 'P3M-AD', 'Reaction-Field'] = 'PME' + coulombtype: Literal["cut-off", "ewald", "PME", "P3M-AD", "Reaction-Field"] = "PME" """ Treatment of electrostatics Allowed values are: @@ -186,7 +191,7 @@ class Config: infinity by setting epsilon-rf =0. Default 'PME' """ - coulomb_modifier: Literal['Potential-shift', 'None'] = 'Potential-shift' + coulomb_modifier: Literal["Potential-shift", "None"] = "Potential-shift" """ Allowed options are: 'Potential-shift': Shift the Coulomb potential by a constant such that it @@ -195,13 +200,13 @@ class Config: comparing energies with those computed with other software. Default 'Potential-shift' """ - rcoulomb_switch: FloatQuantity['nanometer'] = 0 * unit.nanometer + rcoulomb_switch: FloatQuantity["nanometer"] = 0 * unit.nanometer """ Where to start switching the Coulomb potential, only relevant when force or potential switching is used. Default 0 * unit.nanometer """ - rcoulomb: FloatQuantity['nanometer'] = 1.2 * unit.nanometer + rcoulomb: FloatQuantity["nanometer"] = 1.2 * unit.nanometer """" The distance for the Coulomb cut-off. Note that with PME this value can be increased by the PME tuning in gmx mdrun along with the PME grid spacing. @@ -218,7 +223,7 @@ class Config: """ # # # Van der Waals # # # - vdwtype: Literal['Cut-off', 'PME'] = 'Cut-off' + vdwtype: Literal["Cut-off", "PME"] = "Cut-off" """ Treatment of vdW interactions. Allowed options are: 'Cut-off': Plain cut-off with pair list radius rlist and VdW cut-off rvdw, @@ -226,7 +231,9 @@ class Config: 'PME: Fast smooth Particle-mesh Ewald (SPME) for VdW interactions. Default 'Cut-off'. """ - vdw_modifier: Literal['Potential-shift', 'None', 'Force-switch', 'Potential-switch'] = 'Potential-shift' + vdw_modifier: Literal[ + "Potential-shift", "None", "Force-switch", "Potential-switch" + ] = "Potential-shift" """ Allowed values are: 'Potential-shift': Shift the Van der Waals potential by a constant such @@ -244,17 +251,17 @@ class Config: This option should only be used if the force field you are using requires this. """ - rvdw_switch: FloatQuantity['nanometer'] = 0 * unit.nanometer + rvdw_switch: FloatQuantity["nanometer"] = 0 * unit.nanometer """" Where to start switching the LJ force and possibly the potential, only relevant when force or potential switching is used. Default 0 * unit.nanometer """ - rvdw: FloatQuantity['nanometer'] = 1.0 * unit.nanometer + rvdw: FloatQuantity["nanometer"] = 1.0 * unit.nanometer """ Distance for the LJ or Buckingham cut-off. Default 1 * unit.nanometer """ - DispCorr: Literal['no', 'EnerPres', 'Ener'] = 'EnerPres' + DispCorr: Literal["no", "EnerPres", "Ener"] = "EnerPres" """ Allowed values are: 'no': Don’t apply any correction @@ -264,7 +271,7 @@ class Config: """ # # # Ewald # # # - fourierspacing: FloatQuantity['nanometer'] = 0.12 * unit.nanometer + fourierspacing: FloatQuantity["nanometer"] = 0.12 * unit.nanometer """ For ordinary Ewald, the ratio of the box dimensions and the spacing determines a lower bound for the number of wave vectors to use in each @@ -278,9 +285,9 @@ class Config: """ The number of grid points along a dimension to which a charge is mapped. The actual order of the PME interpolation is one less, e.g. the default of - 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for - P3M-AD). When running in parallel, it can be worth to switch to 5 and - simultaneously increase the grid spacing. Note that on the CPU only values + 4 gives cubic interpolation. Supported values are 3 to 12 (max 8 for + P3M-AD). When running in parallel, it can be worth to switch to 5 and + simultaneously increase the grid spacing. Note that on the CPU only values 4 and 5 have SIMD acceleration and GPUs only support the value 4. Default 4. """ @@ -291,7 +298,7 @@ class Config: but then you need more wave vectors for the reciprocal sum. Default 1e-5 """ - lj_pme_comb_rule: Literal['Geometric', 'Lorentz-Berthelot'] = 'Geometric' + lj_pme_comb_rule: Literal["Geometric", "Lorentz-Berthelot"] = "Geometric" """ The combination rules used to combine VdW-parameters in the reciprocal part of LJ-PME. Geometric rules are much faster than Lorentz-Berthelot and @@ -301,13 +308,13 @@ class Config: 'Lorentz-Berthelot': Apply Lorentz-Berthelot combination rules. Default 'Geometric'. """ - ewald_geometry: Literal['3d', '3dc'] = '3d' + ewald_geometry: Literal["3d", "3dc"] = "3d" """ Allowed values are: '3d': The Ewald sum is performed in all three dimensions. '3dc': The reciprocal sum is still performed in 3D, but a force and - potential correction applied in the z dimension to produce a pseudo-2D - summation. + potential correction applied in the z dimension to produce a pseudo-2D + summation. Default '3d'. """ epsilon_surface: float = 0 @@ -321,7 +328,9 @@ class Config: """ # # # Temperature coupling # # # - tcoupl: Literal['no', 'berendsen', 'nose-hoover', 'andersen', 'andersen-massive', 'v-rescale'] = 'no' + tcoupl: Literal[ + "no", "berendsen", "nose-hoover", "andersen", "andersen-massive", "v-rescale" + ] = "no" """ Temperature coupling options. Note that tcoupl will be ignored when the 'sd' integrator is used. @@ -370,23 +379,23 @@ class Config: nsttcouple is set to 1. Default -1. """ - tc_grps: str = 'system' + tc_grps: str = "system" """ Groups to couple to separate temperature baths. Default 'system'. """ - tau_t: FloatQuantity['picosecond'] = 2.0 * unit.picosecond + tau_t: FloatQuantity["picosecond"] = 2.0 * unit.picosecond """" Time constant for coupling (one for each group in tc-grps), -1 means no temperature coupling. Default 2.0 * unit.picosecond. """ - ref_t: FloatQuantity['Kelvin'] = 298.15 * unit.kelvin + ref_t: FloatQuantity["Kelvin"] = 298.15 * unit.kelvin """ Reference temperature for coupling (one for each group in tc-grps). Default 298.15 * unit.kelvin """ # # # Pressure coupling # # # - pcoupl: Literal['no', 'berendsen', 'C-rescale', 'Parrinello-Rahman', 'MTTK'] = 'no' + pcoupl: Literal["no", "berendsen", "C-rescale", "Parrinello-Rahman", "MTTK"] = "no" """ Options for pressure coupling (barostat). Allowed values are: 'no': No pressure coupling. This means a fixed box size. @@ -405,18 +414,20 @@ class Config: As for Nose-Hoover temperature coupling the time constant tau-p is the period of pressure fluctuations at equilibrium. 'MTTK': Martyna-Tuckerman-Tobias-Klein implementation, only useable with - integrator=md-vv or integrator=md-vv-avek, very similar to - Parrinello-Rahman. As for Nose-Hoover temperature coupling the time - constant tau-p is the period of pressure fluctuations at equilibrium. - This is probably a better method when you want to apply pressure - scaling during data collection, but beware that you can get very large - oscillations if you are starting from a different pressure. - This requires a constant ensemble temperature for the system. - Currently it only supports isotropic scaling, and only works without + integrator=md-vv or integrator=md-vv-avek, very similar to + Parrinello-Rahman. As for Nose-Hoover temperature coupling the time + constant tau-p is the period of pressure fluctuations at equilibrium. + This is probably a better method when you want to apply pressure + scaling during data collection, but beware that you can get very large + oscillations if you are starting from a different pressure. + This requires a constant ensemble temperature for the system. + Currently it only supports isotropic scaling, and only works without constraints. Default 'no'. """ - pcoupltype: Literal['isotropic', 'semiisotropic', 'anisotropic', 'surface-tension'] = 'isotropic' + pcoupltype: Literal[ + "isotropic", "semiisotropic", "anisotropic", "surface-tension" + ] = "isotropic" """ Specifies the kind of isotropy of the pressure coupling used. Each kind takes one or more values for compressibility and ref-p. Only a single value @@ -447,30 +458,30 @@ class Config: nstpcouple: int = -1 """ The frequency for coupling the pressure. The default value of -1 sets - nstpcouple equal to 100, or fewer steps if required for accurate - integration (5 steps per tau for first order coupling, 20 steps per tau for - second order coupling). Note that the default value is large in order to - reduce the overhead of the additional computation and communication - required for obtaining the virial and kinetic energy. For velocity Verlet + nstpcouple equal to 100, or fewer steps if required for accurate + integration (5 steps per tau for first order coupling, 20 steps per tau for + second order coupling). Note that the default value is large in order to + reduce the overhead of the additional computation and communication + required for obtaining the virial and kinetic energy. For velocity Verlet integrators nsttcouple is set to 1. Default -1. """ - tau_p: FloatQuantity['picosecond'] = 5 * unit.picosecond + tau_p: FloatQuantity["picosecond"] = 5 * unit.picosecond """ The time constant for pressure coupling (one value for all directions). Default 5 * unit.picosecond. """ - compressibility: FloatQuantity['1/bar'] = 4.5e-05 / unit.bar + compressibility: FloatQuantity["1/bar"] = 4.5e-05 / unit.bar """ The compressibility. For water at 1 atm and 300 K the compressibility is 4.5e-5 bar-1. The number of required values is implied by pcoupltype. Default 4.5e-05 / unit.bar """ - ref_p: FloatQuantity['bar'] = 1.01325 * unit.bar + ref_p: FloatQuantity["bar"] = 1.01325 * unit.bar """ The reference pressure for coupling. The number of required values is implied by pcoupltype. Default 1.01325 * unit.bar. """ - refcoord_scaling: Literal['no', 'all', 'com'] = 'no' + refcoord_scaling: Literal["no", "all", "com"] = "no" """ Allowed values are: 'no': The reference coordinates for position restraints are not modified. @@ -486,17 +497,17 @@ class Config: """ # # # Velocity generation # # # - gen_vel: Literal['no', 'yes'] = 'yes' + gen_vel: Literal["no", "yes"] = "yes" """ Velocity generation. Allowed values are: 'no': Do not generate velocities. The velocities are set to zero when there are no velocities in the input structure file. - 'yes': Generate velocities in gmx grompp according to a Maxwell - distribution at temperature gen-temp, with random seed gen-seed. + 'yes': Generate velocities in gmx grompp according to a Maxwell + distribution at temperature gen-temp, with random seed gen-seed. This is only meaningful with integrator=md. Default 'yes'. """ - gen_temp: FloatQuantity['kelvin'] = 298.15 * unit.kelvin + gen_temp: FloatQuantity["kelvin"] = 298.15 * unit.kelvin """ Temperature for Maxwell distribution. Default 298.15 * unit.kelvin. """ @@ -507,11 +518,13 @@ class Config: """ # # # Bonds # # # - constraints: Literal['none', 'h-bonds', 'all-bonds', 'h-angles', 'all-angles'] = 'h-bonds' + constraints: Literal[ + "none", "h-bonds", "all-bonds", "h-angles", "all-angles" + ] = "h-bonds" """ Controls which bonds in the topology will be converted to rigid holonomic constraints. Note that typical rigid water models do not have bonds, but - rather a specialized [settles] directive, so are not affected by this + rather a specialized [settles] directive, so are not affected by this keyword. Allowed values are: 'none': No bonds converted to constraints. 'h-bonds': Convert the bonds with H-atoms to constraints. @@ -522,7 +535,7 @@ class Config: bond-constraints. Default 'h-bonds' """ - constraint_algorithm: Literal['lincs', 'shake'] = 'lincs' + constraint_algorithm: Literal["lincs", "shake"] = "lincs" """ Chooses which solver satisfies any non-SETTLE holonomic constraints. Allowed values are: @@ -545,7 +558,7 @@ class Config: constrained. SHAKE can not be used with energy minimization. Default 'lincs' """ - continuation: Literal['no', 'yes'] = 'no' + continuation: Literal["no", "yes"] = "no" """ This option was formerly known as unconstrained-start. Allowed values are: @@ -571,12 +584,12 @@ class Config: Number of iterations to correct for rotational lengthening in LINCS. Default 1. """ - lincs_warnangle: FloatQuantity['deg'] = 30 * unit.degree + lincs_warnangle: FloatQuantity["deg"] = 30 * unit.degree """ Maximum angle that a bond can rotate before LINCS will complain. Default 30 * unit.degree. """ - morse: Literal['no', 'yes'] = 'no' + morse: Literal["no", "yes"] = "no" """ Allowed options are: 'no': Bonds are represented by a harmonic potential. @@ -586,9 +599,10 @@ class Config: class OutputSettings(SettingsBaseModel): - """" + """ " Output Settings for simulations run using Gromacs """ + nstxout: int = 0 """ Number of steps that elapse between writing coordinates to the output @@ -658,15 +672,18 @@ class EMSimulationSettings(SimulationSettings): """ Settings for energy minimization. """ - integrator = 'steep' + + integrator = "steep" nsteps = 5000 - emtol: FloatQuantity['kilojoule / (mole * nanometer)'] = 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) + emtol: FloatQuantity["kilojoule / (mole * nanometer)"] = ( + 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) + ) """ The minimization is converged when the maximum force is smaller than this value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) """ - emstep: FloatQuantity['nanometer'] = 0.01 * unit.nanometer + emstep: FloatQuantity["nanometer"] = 0.01 * unit.nanometer """ Initial step size. Default 0.01 * unit.nanometer """ @@ -711,7 +728,7 @@ class Config: Number of independent MD runs to perform. """ - @validator('protocol_repeats') + @validator("protocol_repeats") def must_be_positive(cls, v): if v <= 0: errmsg = f"protocol_repeats must be a positive value, got {v}." From 8c8c95e3404e952686fdb1833f261ee7193864a7 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Tue, 16 Jul 2024 12:53:54 +0200 Subject: [PATCH 08/23] small fix --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 65c8372..c3d362c 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -7,7 +7,7 @@ :class:`openfe.protocols.gromacs_md.md_methods.py` """ -from typing import Literal, Optional +from typing import Literal from gufe.settings import OpenMMSystemGeneratorFFSettings, SettingsBaseModel from openfe.protocols.openmm_utils.omm_settings import ( From 89486dcf5356acf7c5057f971ce5dc0a825a5ec3 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Tue, 16 Jul 2024 13:22:25 +0200 Subject: [PATCH 09/23] Add details EM, NVT, NPT --- .../protocols/gromacs_md/md_settings.py | 22 +++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index c3d362c..1d702e3 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -7,7 +7,7 @@ :class:`openfe.protocols.gromacs_md.md_methods.py` """ -from typing import Literal +from typing import Literal, Optional from gufe.settings import OpenMMSystemGeneratorFFSettings, SettingsBaseModel from openfe.protocols.openmm_utils.omm_settings import ( @@ -103,7 +103,7 @@ class Config: Frequency for center of mass motion removal in unit [steps]. Default 100 """ - comm_grps: str + comm_grps: Optional[str] """ Group(s) for center of mass motion removal, default is the whole system. """ @@ -656,12 +656,12 @@ class OutputSettings(SettingsBaseModel): Precision with which to write to the compressed trajectory file. Default 1000. """ - compressed_x_grps: str + compressed_x_grps: Optional[str] """ Group(s) to write to the compressed trajectory file, by default the whole system is written (if nstxout-compressed > 0). """ - energygrps: str + energygrps: Optional[str] """ Group(s) for which to write to write short-ranged non-bonded potential energies to the energy file (not supported on GPUs) @@ -687,6 +687,9 @@ class EMSimulationSettings(SimulationSettings): """ Initial step size. Default 0.01 * unit.nanometer """ + tcoupl = 'no' + pcoupl = 'no' + gen_vel = 'no' class EMOutputSettings(OutputSettings): @@ -699,6 +702,9 @@ class NVTSimulationSettings(SimulationSettings): """ Settings for MD simulation in the NVT ensemble. """ + nsteps = 50000 # 100ps + pcoupl = 'no' + gen_vel = 'yes' class NVTOutputSettings(OutputSettings): @@ -711,12 +717,20 @@ class NPTSimulationSettings(SimulationSettings): """ Settings for MD simulation in the NPT ensemble. """ + pcoupl = 'Parrinello-Rahman' + pcoupltype = 'isotropic' + ref_p = 1.01325 * unit.bar + gen_vel = 'no' # If continuation from NVT simulation class NPTOutputSettings(OutputSettings): """ Output Settings for the MD simulation in the NPT ensemble. """ + nstxout = 5000 + nstvout = 5000 + nstfout = 5000 + nstxout_compressed = 5000 class GromacsMDProtocolSettings(Settings): From cbedd1a591ef5d01a1cc7e8fb4f700a363b57fee Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 16 Jul 2024 11:22:38 +0000 Subject: [PATCH 10/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../protocols/gromacs_md/md_settings.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 1d702e3..4b37c2e 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -687,9 +687,9 @@ class EMSimulationSettings(SimulationSettings): """ Initial step size. Default 0.01 * unit.nanometer """ - tcoupl = 'no' - pcoupl = 'no' - gen_vel = 'no' + tcoupl = "no" + pcoupl = "no" + gen_vel = "no" class EMOutputSettings(OutputSettings): @@ -702,9 +702,10 @@ class NVTSimulationSettings(SimulationSettings): """ Settings for MD simulation in the NVT ensemble. """ + nsteps = 50000 # 100ps - pcoupl = 'no' - gen_vel = 'yes' + pcoupl = "no" + gen_vel = "yes" class NVTOutputSettings(OutputSettings): @@ -717,16 +718,18 @@ class NPTSimulationSettings(SimulationSettings): """ Settings for MD simulation in the NPT ensemble. """ - pcoupl = 'Parrinello-Rahman' - pcoupltype = 'isotropic' + + pcoupl = "Parrinello-Rahman" + pcoupltype = "isotropic" ref_p = 1.01325 * unit.bar - gen_vel = 'no' # If continuation from NVT simulation + gen_vel = "no" # If continuation from NVT simulation class NPTOutputSettings(OutputSettings): """ Output Settings for the MD simulation in the NPT ensemble. """ + nstxout = 5000 nstvout = 5000 nstfout = 5000 From ac39a257bf1719287373d0b147732a28b1f2de84 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Tue, 16 Jul 2024 15:06:57 +0200 Subject: [PATCH 11/23] small fixes2 --- .../protocols/gromacs_md/md_settings.py | 29 ++++++------------- 1 file changed, 9 insertions(+), 20 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 4b37c2e..635bb79 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -603,6 +603,11 @@ class OutputSettings(SettingsBaseModel): Output Settings for simulations run using Gromacs """ + forcefield_cache: Optional[str] = 'db.json' + """ + Filename for caching small molecule residue templates so they can be + later reused. + """ nstxout: int = 0 """ Number of steps that elapse between writing coordinates to the output @@ -673,9 +678,6 @@ class EMSimulationSettings(SimulationSettings): Settings for energy minimization. """ - integrator = "steep" - nsteps = 5000 - emtol: FloatQuantity["kilojoule / (mole * nanometer)"] = ( 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) ) @@ -687,9 +689,6 @@ class EMSimulationSettings(SimulationSettings): """ Initial step size. Default 0.01 * unit.nanometer """ - tcoupl = "no" - pcoupl = "no" - gen_vel = "no" class EMOutputSettings(OutputSettings): @@ -703,10 +702,6 @@ class NVTSimulationSettings(SimulationSettings): Settings for MD simulation in the NVT ensemble. """ - nsteps = 50000 # 100ps - pcoupl = "no" - gen_vel = "yes" - class NVTOutputSettings(OutputSettings): """ @@ -719,22 +714,12 @@ class NPTSimulationSettings(SimulationSettings): Settings for MD simulation in the NPT ensemble. """ - pcoupl = "Parrinello-Rahman" - pcoupltype = "isotropic" - ref_p = 1.01325 * unit.bar - gen_vel = "no" # If continuation from NVT simulation - class NPTOutputSettings(OutputSettings): """ Output Settings for the MD simulation in the NPT ensemble. """ - nstxout = 5000 - nstvout = 5000 - nstfout = 5000 - nstxout_compressed = 5000 - class GromacsMDProtocolSettings(Settings): class Config: @@ -752,6 +737,10 @@ def must_be_positive(cls, v): raise ValueError(errmsg) return v + # File names for .gro and .top file + gro: str + top: str + # Things for creating the systems forcefield_settings: OpenMMSystemGeneratorFFSettings partial_charge_settings: OpenFFPartialChargeSettings From 5824192ae93c2c003d04b124748e6dfa62422d0e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 16 Jul 2024 13:07:35 +0000 Subject: [PATCH 12/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 635bb79..4600ab5 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -603,7 +603,7 @@ class OutputSettings(SettingsBaseModel): Output Settings for simulations run using Gromacs """ - forcefield_cache: Optional[str] = 'db.json' + forcefield_cache: Optional[str] = "db.json" """ Filename for caching small molecule residue templates so they can be later reused. From 413ce394f5dea546a90e64aa5e36e44b24636944 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Tue, 23 Jul 2024 13:03:35 +0200 Subject: [PATCH 13/23] Add initial validators --- .../protocols/gromacs_md/md_settings.py | 67 +++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 4600ab5..8326159 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -597,6 +597,64 @@ class Config: Default 'no'. """ + @validator( + "tinit", + "nsteps", + "init_step", + "simulation_part", + "nstcomm", + "nstlist", + "rlist", + "rcoulomb_switch", + "rcoulomb", + "epsilon_r", + "epsilon_rf", + "rvdw_switch", + "rvdw", + "fourierspacing", + "ewald_rtol", + "epsilon_surface", + "tau_t", + "ref_t", + "tau_p", + "compressibility", + "gen_temp", + "shake_tol", + "lincs_order", + "lincs_iter", + "lincs_warnangle", + ) + def must_be_positive_or_zero(cls, v): + if v < 0: + errmsg = ( + "Settings tinit, nsteps, init_step, simulation_part, " + "nstcomm, nstlist, rlist, rcoulomb_switch, rcoulomb, " + "epsilon_r, epsilon_rf, rvdw_switch, rvdw, " + "fourierspacing, ewald_rtol, epsilon_surface, tau_t, " + "ref_t, tau_p, compressibility, gen_temp, shake_tol, " + "lincs_order, lincs_iter, and lincs_warnangle must be" + f" zero or positive values, got {v}." + ) + raise ValueError(errmsg) + return v + + @validator("dt", "mass_repartition_factor") + def must_be_positive(cls, v): + if v <= 0: + errmsg = ( + "timestep dt, and mass_repartition_factor " + f"must be positive values, got {v}." + ) + raise ValueError(errmsg) + return v + + @validator("pme_order") + def must_be_between_3_12(cls, v): + if not 3 <= v <= 12: + errmsg = "pme_order " f"must be between 3 and 12, got {v}." + raise ValueError(errmsg) + return v + class OutputSettings(SettingsBaseModel): """ " @@ -690,6 +748,15 @@ class EMSimulationSettings(SimulationSettings): Initial step size. Default 0.01 * unit.nanometer """ + @validator('integrator') + def is_steep(cls, v): + # EM should have 'steep' integrator + if v != 'steep': + errmsg = ("For energy minimization, only the integrator=steep " + f"is supported, got integrator={v}.") + raise ValueError(errmsg) + return v + class EMOutputSettings(OutputSettings): """ From 37f1efa19b9bfbcef21e8448db1f4d70a28f1d45 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 23 Jul 2024 11:04:01 +0000 Subject: [PATCH 14/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 8326159..ca51121 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -748,15 +748,17 @@ class EMSimulationSettings(SimulationSettings): Initial step size. Default 0.01 * unit.nanometer """ - @validator('integrator') + @validator("integrator") def is_steep(cls, v): # EM should have 'steep' integrator - if v != 'steep': - errmsg = ("For energy minimization, only the integrator=steep " - f"is supported, got integrator={v}.") + if v != "steep": + errmsg = ( + "For energy minimization, only the integrator=steep " + f"is supported, got integrator={v}." + ) raise ValueError(errmsg) return v - + class EMOutputSettings(OutputSettings): """ From 8196443d019f849112a6b438d17635b37d4d8791 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Wed, 24 Jul 2024 15:07:31 +0200 Subject: [PATCH 15/23] Update Settings with changes from other PR --- .../protocols/gromacs_md/md_settings.py | 390 +++++------------- 1 file changed, 111 insertions(+), 279 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index ca51121..89d93d7 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -35,78 +35,31 @@ class Config: arbitrary_types_allowed = True # # # Run control # # # - integrator: Literal["md", "md-vv", "md-vv-avek", "sd", "steep"] = "sd" + integrator: Literal["md", "sd", "steep"] = "sd" """ MD integrators and other algorithms (steep for energy minimization). Allowed values are: 'md': a leap-frog integrator - 'md-vv': a velocity Verlet integrator - 'md-vv-avek': A velocity Verlet algorithm identical to integrator=md-vv, - except that the kinetic energy is determined as the average - of the two half step kinetic energies as in the - integrator=md integrator, and this thus more accurate. 'sd': A leap-frog stochastic dynamics integrator. Note that when using this integrator the parameters tcoupl and nsttcouple are ignored. 'steep': A steepest descent algorithm for energy minimization. """ - tinit: FloatQuantity["picosecond"] = 0 * unit.picosecond - """ - Starting time for the MD run. - This only makes sense for time-based integrators. - Default 0 * unit.picosecond - """ - dt: FloatQuantity["picosecond"] = 0.001 * unit.picosecond + dt: FloatQuantity["picosecond"] = 0.002 * unit.picosecond """ Time step for integration (only makes sense for time-based integrators). - Default 0.001 * unit.picosecond + Default 0.002 * unit.picosecond """ nsteps: int = 0 """ Maximum number of steps to integrate or minimize, -1 is no maximum Default 0 """ - init_step: int = 0 - """ - The starting step. The time at step i in a run is calculated as: - t = tinit + dt * (init-step + i). - Default 0 - """ - simulation_part: int = 0 - """ - A simulation can consist of multiple parts, each of which has a part number. - This option specifies what that number will be, which helps keep track of - parts that are logically the same simulation. This option is generally - useful to set only when coping with a crashed simulation where files were - lost. Default 0 - """ mass_repartition_factor: int = 1 """ Scales the masses of the lightest atoms in the system by this factor to the mass mMin. All atoms with a mass lower than mMin also have their mass set to that mMin. Default 1 (no mass scaling) """ - comm_mode: Literal[ - "Linear", "Angular", "Linear-acceleration-correction", "None" - ] = "Linear" - """ - Settings for center of mass treatmeant - Allowed values are: - 'Linear': Remove center of mass translational velocity - 'Angular': Remove center of mass translational and rotational velocity - 'Linear-acceleration-correction': Remove center of mass translational - velocity. Correct the center of mass position assuming linear acceleration - over nstcomm steps. - 'None': No restriction on the center of mass motion - """ - nstcomm: int = 100 - """ - Frequency for center of mass motion removal in unit [steps]. - Default 100 - """ - comm_grps: Optional[str] - """ - Group(s) for center of mass motion removal, default is the whole system. - """ # # # Langevin dynamics # # # ld_seed: int = -1 @@ -142,24 +95,6 @@ class Config: This can be used in combination with walls. Default 'xyz'. """ - verlet_buffer_tolerance: FloatQuantity["kilojoule / (mole * picosecond)"] = ( - 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) - ) - """" - Used when performing a simulation with dynamics. This sets the maximum - allowed error for pair interactions per particle caused by the Verlet - buffer, which indirectly sets rlist. - Default 0.005 * unit.kilojoule / (unit.mole * unit.picosecond) - """ - verlet_buffer_pressure_tolerance: FloatQuantity["bar"] = 0.5 * unit.bar - """ - Used when performing a simulation with dynamics and only active when - verlet-buffer-tolerance is positive. This sets the maximum tolerated error - in the average pressure due to missing Lennard-Jones interactions of - particle pairs that are not in the pair list, but come within rvdw range - as the pair list ages. - Default 0.5 * unit.bar - """ rlist: FloatQuantity["nanometer"] = 1 * unit.nanometer """ Cut-off distance for the short-range neighbor list. With dynamics, this is @@ -191,36 +126,12 @@ class Config: infinity by setting epsilon-rf =0. Default 'PME' """ - coulomb_modifier: Literal["Potential-shift", "None"] = "Potential-shift" - """ - Allowed options are: - 'Potential-shift': Shift the Coulomb potential by a constant such that it - is zero at the cut-off. - 'None': Use an unmodified Coulomb potential. This can be useful when - comparing energies with those computed with other software. - Default 'Potential-shift' - """ - rcoulomb_switch: FloatQuantity["nanometer"] = 0 * unit.nanometer - """ - Where to start switching the Coulomb potential, only relevant when force - or potential switching is used. - Default 0 * unit.nanometer - """ rcoulomb: FloatQuantity["nanometer"] = 1.2 * unit.nanometer """" The distance for the Coulomb cut-off. Note that with PME this value can be increased by the PME tuning in gmx mdrun along with the PME grid spacing. Default 1.2 * unit.nanometer """ - epsilon_r: int = 1 - """ - The relative dielectric constant. A value of 0 means infinity. Default 1. - """ - epsilon_rf: int = 0 - """ - The relative dielectric constant of the reaction field. This is only used - with reaction-field electrostatics. A value of 0 means infinity. Default 0 - """ # # # Van der Waals # # # vdwtype: Literal["Cut-off", "PME"] = "Cut-off" @@ -231,31 +142,13 @@ class Config: 'PME: Fast smooth Particle-mesh Ewald (SPME) for VdW interactions. Default 'Cut-off'. """ - vdw_modifier: Literal[ - "Potential-shift", "None", "Force-switch", "Potential-switch" - ] = "Potential-shift" + vdw_modifier: Literal["Potential-shift", "None"] = "Potential-shift" """ Allowed values are: 'Potential-shift': Shift the Van der Waals potential by a constant such that it is zero at the cut-off. 'None': Use an unmodified Van der Waals potential. This can be useful when comparing energies with those computed with other software. - 'Force-switch': Smoothly switches the forces to zero between rvdw-switch - and rvdw. This shifts the potential shift over the whole range and - switches it to zero at the cut-off. Note that this is more expensive - to calculate than a plain cut-off and it is not required for energy - conservation, since Potential-shift conserves energy just as well. - 'Potential-switch': Smoothly switches the potential to zero between - rvdw-switch and rvdw. Note that this introduces articifically large - forces in the switching region and is much more expensive to calculate. - This option should only be used if the force field you are using - requires this. - """ - rvdw_switch: FloatQuantity["nanometer"] = 0 * unit.nanometer - """" - Where to start switching the LJ force and possibly the potential, only - relevant when force or potential switching is used. - Default 0 * unit.nanometer """ rvdw: FloatQuantity["nanometer"] = 1.0 * unit.nanometer """ @@ -271,16 +164,6 @@ class Config: """ # # # Ewald # # # - fourierspacing: FloatQuantity["nanometer"] = 0.12 * unit.nanometer - """ - For ordinary Ewald, the ratio of the box dimensions and the spacing - determines a lower bound for the number of wave vectors to use in each - (signed) direction. For PME and P3M, that ratio determines a lower bound - for the number of Fourier-space grid points that will be used along that - axis. In all cases, the number for each direction can be overridden by - entering a non-zero value for that fourier-nx direction. - Default 0.12 * unit.nanometer - """ pme_order: int = 4 """ The number of grid points along a dimension to which a charge is mapped. @@ -298,34 +181,6 @@ class Config: but then you need more wave vectors for the reciprocal sum. Default 1e-5 """ - lj_pme_comb_rule: Literal["Geometric", "Lorentz-Berthelot"] = "Geometric" - """ - The combination rules used to combine VdW-parameters in the reciprocal part - of LJ-PME. Geometric rules are much faster than Lorentz-Berthelot and - usually the recommended choice, even when the rest of the force field uses - the Lorentz-Berthelot rules. Allowed values are: - 'Geometric': Apply geometric combination rules. - 'Lorentz-Berthelot': Apply Lorentz-Berthelot combination rules. - Default 'Geometric'. - """ - ewald_geometry: Literal["3d", "3dc"] = "3d" - """ - Allowed values are: - '3d': The Ewald sum is performed in all three dimensions. - '3dc': The reciprocal sum is still performed in 3D, but a force and - potential correction applied in the z dimension to produce a pseudo-2D - summation. - Default '3d'. - """ - epsilon_surface: float = 0 - """ - This controls the dipole correction to the Ewald summation in 3D. - The default value of zero means it is turned off. Turn it on by setting it - to the value of the relative permittivity of the imaginary surface around - your infinite system. Be careful - you shouldn’t use this if you have free - mobile charges in your system. This value does not affect the slab 3DC - variant of the long range corrections. Default 0. - """ # # # Temperature coupling # # # tcoupl: Literal[ @@ -368,34 +223,14 @@ class Config: is written to the energy and log file. Default 'no' (for the default integrator, 'sd', this option is ignored). """ - nsttcouple: int = -1 - """ - The frequency for coupling the temperature. The default value of -1 sets - nsttcouple equal to 100, or fewer steps if required for accurate - integration (5 steps per tau for first order coupling, 20 steps per tau for - second order coupling). Note that the default value is large in order to - reduce the overhead of the additional computation and communication - required for obtaining the kinetic energy. For velocity Verlet integrators - nsttcouple is set to 1. - Default -1. - """ - tc_grps: str = "system" - """ - Groups to couple to separate temperature baths. Default 'system'. - """ - tau_t: FloatQuantity["picosecond"] = 2.0 * unit.picosecond - """" - Time constant for coupling (one for each group in tc-grps), -1 means no - temperature coupling. Default 2.0 * unit.picosecond. - """ - ref_t: FloatQuantity["Kelvin"] = 298.15 * unit.kelvin + ref_t: FloatQuantity["kelvin"] = 298.15 * unit.kelvin """ Reference temperature for coupling (one for each group in tc-grps). Default 298.15 * unit.kelvin """ # # # Pressure coupling # # # - pcoupl: Literal["no", "berendsen", "C-rescale", "Parrinello-Rahman", "MTTK"] = "no" + pcoupl: Literal["no", "berendsen", "C-rescale", "Parrinello-Rahman"] = "no" """ Options for pressure coupling (barostat). Allowed values are: 'no': No pressure coupling. This means a fixed box size. @@ -413,69 +248,8 @@ class Config: for the atoms is coupled to this. No instantaneous scaling takes place. As for Nose-Hoover temperature coupling the time constant tau-p is the period of pressure fluctuations at equilibrium. - 'MTTK': Martyna-Tuckerman-Tobias-Klein implementation, only useable with - integrator=md-vv or integrator=md-vv-avek, very similar to - Parrinello-Rahman. As for Nose-Hoover temperature coupling the time - constant tau-p is the period of pressure fluctuations at equilibrium. - This is probably a better method when you want to apply pressure - scaling during data collection, but beware that you can get very large - oscillations if you are starting from a different pressure. - This requires a constant ensemble temperature for the system. - Currently it only supports isotropic scaling, and only works without - constraints. Default 'no'. """ - pcoupltype: Literal[ - "isotropic", "semiisotropic", "anisotropic", "surface-tension" - ] = "isotropic" - """ - Specifies the kind of isotropy of the pressure coupling used. Each kind - takes one or more values for compressibility and ref-p. Only a single value - is permitted for tau-p. Allowed values are: - 'isotropic': sotropic pressure coupling with time constant tau-p. One value - each for compressibility and ref-p is required. - 'semiisotropic': Pressure coupling which is isotropic in the x and y - direction, but different in the z direction. This can be useful for - membrane simulations. Two values each for compressibility and ref-p are - required, for x/y and z directions respectively. - 'anisotropic': Same as before, but 6 values are needed for xx, yy, zz, - xy/yx, xz/zx and yz/zy components, respectively. When the off-diagonal - compressibilities are set to zero, a rectangular box will stay - rectangular. Beware that anisotropic scaling can lead to extreme - deformation of the simulation box. - 'surface-tension': Surface tension coupling for surfaces parallel to the - xy-plane. Uses normal pressure coupling for the z-direction, while the - surface tension is coupled to the x/y dimensions of the box. The first - ref-p value is the reference surface tension times the number of - surfaces bar nm, the second value is the reference z-pressure bar. - The two compressibility values are the compressibility in the x/y and z - direction respectively. The value for the z-compressibility should be - reasonably accurate since it influences the convergence of the - surface-tension, it can also be set to zero to have a box with constant - height. - Default 'isotropic' - """ - nstpcouple: int = -1 - """ - The frequency for coupling the pressure. The default value of -1 sets - nstpcouple equal to 100, or fewer steps if required for accurate - integration (5 steps per tau for first order coupling, 20 steps per tau for - second order coupling). Note that the default value is large in order to - reduce the overhead of the additional computation and communication - required for obtaining the virial and kinetic energy. For velocity Verlet - integrators nsttcouple is set to 1. Default -1. - """ - tau_p: FloatQuantity["picosecond"] = 5 * unit.picosecond - """ - The time constant for pressure coupling (one value for all directions). - Default 5 * unit.picosecond. - """ - compressibility: FloatQuantity["1/bar"] = 4.5e-05 / unit.bar - """ - The compressibility. For water at 1 atm and 300 K the compressibility is - 4.5e-5 bar-1. The number of required values is implied by pcoupltype. - Default 4.5e-05 / unit.bar - """ ref_p: FloatQuantity["bar"] = 1.01325 * unit.bar """ The reference pressure for coupling. The number of required values is @@ -584,56 +358,26 @@ class Config: Number of iterations to correct for rotational lengthening in LINCS. Default 1. """ - lincs_warnangle: FloatQuantity["deg"] = 30 * unit.degree - """ - Maximum angle that a bond can rotate before LINCS will complain. - Default 30 * unit.degree. - """ - morse: Literal["no", "yes"] = "no" - """ - Allowed options are: - 'no': Bonds are represented by a harmonic potential. - 'yes': Bonds are represented by a Morse potential. - Default 'no'. - """ @validator( - "tinit", "nsteps", - "init_step", - "simulation_part", - "nstcomm", "nstlist", "rlist", - "rcoulomb_switch", "rcoulomb", - "epsilon_r", - "epsilon_rf", - "rvdw_switch", "rvdw", - "fourierspacing", "ewald_rtol", - "epsilon_surface", - "tau_t", "ref_t", - "tau_p", - "compressibility", "gen_temp", "shake_tol", "lincs_order", "lincs_iter", - "lincs_warnangle", ) def must_be_positive_or_zero(cls, v): if v < 0: errmsg = ( - "Settings tinit, nsteps, init_step, simulation_part, " - "nstcomm, nstlist, rlist, rcoulomb_switch, rcoulomb, " - "epsilon_r, epsilon_rf, rvdw_switch, rvdw, " - "fourierspacing, ewald_rtol, epsilon_surface, tau_t, " - "ref_t, tau_p, compressibility, gen_temp, shake_tol, " - "lincs_order, lincs_iter, and lincs_warnangle must be" - f" zero or positive values, got {v}." + "Settings nsteps, nstlist, rlist, rcoulomb, rvdw, ewald_rtol, " + "ref_t, gen_temp, shake_tol, lincs_order, and lincs_iter" + f" must be zero or positive values, got {v}." ) raise ValueError(errmsg) return v @@ -655,6 +399,42 @@ def must_be_between_3_12(cls, v): raise ValueError(errmsg) return v + @validator("dt") + def is_time(cls, v): + if not v.is_compatible_with(unit.picosecond): + raise ValueError("dt must be in time units " "(i.e. picoseconds)") + + @validator("rlist", "rcoulomb", "rvdw") + def is_distance(cls, v): + if not v.is_compatible_with(unit.nanometer): + raise ValueError( + "rlist, rcoulomb, and rvdw must be in distance " + "units (i.e. nanometers)" + ) + + @validator("ref_t", "gen_temp") + def is_temperature(cls, v): + if not v.is_compatible_with(unit.kelvin): + raise ValueError( + "ref_t and gen_temp must be in temperature units" " (i.e. kelvin)" + ) + + @validator("ref_p") + def is_pressure(cls, v): + if not v.is_compatible_with(unit.bar): + raise ValueError("ref_p must be in pressure units (i.e. bar)") + + @validator("integrator") + def supported_integrator(cls, v): + supported = ["md", "sd", "steep"] + if v.lower() not in supported: + errmsg = ( + "Only the following sampler_method values are " + f"supported: {supported}, got {v}" + ) + raise ValueError(errmsg) + return v + class OutputSettings(SettingsBaseModel): """ " @@ -666,6 +446,11 @@ class OutputSettings(SettingsBaseModel): Filename for caching small molecule residue templates so they can be later reused. """ + mdp_file: str = "em.mdp" + """ + Filename for the mdp file for running simulations in Gromacs. + Default 'em.mdp' + """ nstxout: int = 0 """ Number of steps that elapse between writing coordinates to the output @@ -719,35 +504,44 @@ class OutputSettings(SettingsBaseModel): Precision with which to write to the compressed trajectory file. Default 1000. """ - compressed_x_grps: Optional[str] + compressed_x_grps: str = "" """ Group(s) to write to the compressed trajectory file, by default the whole system is written (if nstxout-compressed > 0). """ - energygrps: Optional[str] + energygrps: str = "" """ Group(s) for which to write to write short-ranged non-bonded potential energies to the energy file (not supported on GPUs) """ + @validator( + "nstxout", + "nstvout", + "nstfout", + "nstlog", + "nstcalcenergy", + "nstenergy", + "nstxout_compressed", + "compressed_x_precision", + ) + def must_be_positive_or_zero(cls, v): + if v < 0: + errmsg = ( + "nstxout, nstvout, nstfout, nstlog, nstcalcenergy, " + "nstenergy, nstxout_compressed, and " + "compressed_x_precision must be zero or positive values," + f" got {v}." + ) + raise ValueError(errmsg) + return v + class EMSimulationSettings(SimulationSettings): """ Settings for energy minimization. """ - emtol: FloatQuantity["kilojoule / (mole * nanometer)"] = ( - 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) - ) - """ - The minimization is converged when the maximum force is smaller than this - value. Default 10.0 * unit.kilojoule / (unit.mole * unit.nanometer) - """ - emstep: FloatQuantity["nanometer"] = 0.01 * unit.nanometer - """ - Initial step size. Default 0.01 * unit.nanometer - """ - @validator("integrator") def is_steep(cls, v): # EM should have 'steep' integrator @@ -771,6 +565,25 @@ class NVTSimulationSettings(SimulationSettings): Settings for MD simulation in the NVT ensemble. """ + @validator("integrator") + def is_not_steep(cls, v): + # needs an MD integrator + if v == "steep": + errmsg = ( + "Molecular Dynamics settings need an MD integrator, " + f"not integrator={v}." + ) + raise ValueError(errmsg) + return v + + @validator("pcoupl") + def has_no_barostat(cls, v): + # NVT cannot have a barostat + if v != "no": + errmsg = f"NVT settings cannot have a barostat, got pcoupl={v}." + raise ValueError(errmsg) + return v + class NVTOutputSettings(OutputSettings): """ @@ -783,6 +596,25 @@ class NPTSimulationSettings(SimulationSettings): Settings for MD simulation in the NPT ensemble. """ + @validator("integrator") + def is_not_steep(cls, v): + # needs an MD integrator + if v == "steep": + errmsg = ( + "Molecular Dynamics settings need an MD integrator, " + f"not integrator={v}." + ) + raise ValueError(errmsg) + return v + + @validator("pcoupl") + def has_barostat(cls, v): + # NPT needs a barostat + if v == "no": + errmsg = f"NPT settings need a barostat, got pcoupl={v}." + raise ValueError(errmsg) + return v + class NPTOutputSettings(OutputSettings): """ From c462d6162ff62b321badd91da1955111e6dedf20 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Thu, 1 Aug 2024 14:03:48 +0200 Subject: [PATCH 16/23] Add returns for validators --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 89d93d7..853f4d7 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -403,6 +403,7 @@ def must_be_between_3_12(cls, v): def is_time(cls, v): if not v.is_compatible_with(unit.picosecond): raise ValueError("dt must be in time units " "(i.e. picoseconds)") + return v @validator("rlist", "rcoulomb", "rvdw") def is_distance(cls, v): @@ -411,6 +412,7 @@ def is_distance(cls, v): "rlist, rcoulomb, and rvdw must be in distance " "units (i.e. nanometers)" ) + return v @validator("ref_t", "gen_temp") def is_temperature(cls, v): @@ -418,11 +420,13 @@ def is_temperature(cls, v): raise ValueError( "ref_t and gen_temp must be in temperature units" " (i.e. kelvin)" ) + return v @validator("ref_p") def is_pressure(cls, v): if not v.is_compatible_with(unit.bar): raise ValueError("ref_p must be in pressure units (i.e. bar)") + return v @validator("integrator") def supported_integrator(cls, v): From 9eb492228b246b74a08e893f8ccd5d17297374ae Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Fri, 2 Aug 2024 11:33:00 +0200 Subject: [PATCH 17/23] Small fix --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 853f4d7..904803f 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -4,7 +4,7 @@ """Settings class for plain MD Protocols using Gromacs + OpenMMTools This module implements the settings necessary to run MD simulations using -:class:`openfe.protocols.gromacs_md.md_methods.py` +:class:`openfe_gromacs.protocols.gromacs_md.md_methods.py` """ from typing import Literal, Optional From 6d65508783ead795bb081b94d6987b32821a9981 Mon Sep 17 00:00:00 2001 From: Hannah Baumann <43765638+hannahbaumann@users.noreply.github.com> Date: Wed, 4 Sep 2024 13:37:10 +0200 Subject: [PATCH 18/23] Update openfe_gromacs/protocols/gromacs_md/md_settings.py Co-authored-by: Irfan Alibay --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 904803f..ac08f33 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -20,10 +20,7 @@ from openff.models.types import FloatQuantity from openff.units import unit -try: - from pydantic.v1 import validator -except ImportError: - from pydantic import validator # type: ignore[assignment] +from pydantic.v1 import validator class SimulationSettings(SettingsBaseModel): From 876eab9b220739c69fc9587c55e275b081e8bc10 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 4 Sep 2024 11:37:21 +0000 Subject: [PATCH 19/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- openfe_gromacs/protocols/gromacs_md/md_settings.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index ac08f33..91cb0ae 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -19,7 +19,6 @@ ) from openff.models.types import FloatQuantity from openff.units import unit - from pydantic.v1 import validator From e141ad80bd2ce596ac224b73c4ae989a64f62d5a Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Wed, 4 Sep 2024 15:58:29 +0200 Subject: [PATCH 20/23] Move integrator to EM and MD Settings --- .../protocols/gromacs_md/md_settings.py | 76 +++++++------------ 1 file changed, 28 insertions(+), 48 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 91cb0ae..985b24a 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -31,20 +31,6 @@ class Config: arbitrary_types_allowed = True # # # Run control # # # - integrator: Literal["md", "sd", "steep"] = "sd" - """ - MD integrators and other algorithms (steep for energy minimization). - Allowed values are: - 'md': a leap-frog integrator - 'sd': A leap-frog stochastic dynamics integrator. Note that when using this - integrator the parameters tcoupl and nsttcouple are ignored. - 'steep': A steepest descent algorithm for energy minimization. - """ - dt: FloatQuantity["picosecond"] = 0.002 * unit.picosecond - """ - Time step for integration (only makes sense for time-based integrators). - Default 0.002 * unit.picosecond - """ nsteps: int = 0 """ Maximum number of steps to integrate or minimize, -1 is no maximum @@ -542,16 +528,12 @@ class EMSimulationSettings(SimulationSettings): Settings for energy minimization. """ - @validator("integrator") - def is_steep(cls, v): - # EM should have 'steep' integrator - if v != "steep": - errmsg = ( - "For energy minimization, only the integrator=steep " - f"is supported, got integrator={v}." - ) - raise ValueError(errmsg) - return v + integrator: Literal["steep"] = "steep" + """ + Method for energy minimization. + Allowed value is: + 'steep': A steepest descent algorithm for energy minimization. + """ class EMOutputSettings(OutputSettings): @@ -560,21 +542,30 @@ class EMOutputSettings(OutputSettings): """ -class NVTSimulationSettings(SimulationSettings): +class MDSimulationSettings(SimulationSettings): """ - Settings for MD simulation in the NVT ensemble. + Settings for MD simulations """ - @validator("integrator") - def is_not_steep(cls, v): - # needs an MD integrator - if v == "steep": - errmsg = ( - "Molecular Dynamics settings need an MD integrator, " - f"not integrator={v}." - ) - raise ValueError(errmsg) - return v + integrator: Literal["md", "sd"] = "sd" + """ + MD integrators and other algorithms (steep for energy minimization). + Allowed values are: + 'md': a leap-frog integrator + 'sd': A leap-frog stochastic dynamics integrator. Note that when using this + integrator the parameters tcoupl and nsttcouple are ignored. + """ + dt: FloatQuantity["picosecond"] = 0.002 * unit.picosecond + """ + Time step for integration (only makes sense for time-based integrators). + Default 0.002 * unit.picosecond + """ + + +class NVTSimulationSettings(MDSimulationSettings): + """ + Settings for MD simulation in the NVT ensemble. + """ @validator("pcoupl") def has_no_barostat(cls, v): @@ -591,22 +582,11 @@ class NVTOutputSettings(OutputSettings): """ -class NPTSimulationSettings(SimulationSettings): +class NPTSimulationSettings(MDSimulationSettings): """ Settings for MD simulation in the NPT ensemble. """ - @validator("integrator") - def is_not_steep(cls, v): - # needs an MD integrator - if v == "steep": - errmsg = ( - "Molecular Dynamics settings need an MD integrator, " - f"not integrator={v}." - ) - raise ValueError(errmsg) - return v - @validator("pcoupl") def has_barostat(cls, v): # NPT needs a barostat From 52aa17ed1faf31a9bd70ff4fd42604f4513b3fd3 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Wed, 4 Sep 2024 16:37:23 +0200 Subject: [PATCH 21/23] Addressing review comments --- .../protocols/gromacs_md/md_settings.py | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 985b24a..43da700 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -67,16 +67,6 @@ class Config: If set to zero the neighbor list is only constructed once and never updated. Default 10. """ - pbc: Literal["xyz", "no", "xy"] = "xyz" - """ - Treatment of periodic boundary conditions. - Allowed values are: - 'xyz': Use periodic boundary conditions in all directions. - 'no': Use no periodic boundary conditions, ignore the box. - 'xy': Use periodic boundary conditions in x and y directions only. - This can be used in combination with walls. - Default 'xyz'. - """ rlist: FloatQuantity["nanometer"] = 1 * unit.nanometer """ Cut-off distance for the short-range neighbor list. With dynamics, this is @@ -314,14 +304,6 @@ class Config: constrained. SHAKE can not be used with energy minimization. Default 'lincs' """ - continuation: Literal["no", "yes"] = "no" - """ - This option was formerly known as unconstrained-start. - Allowed values are: - 'no': Apply constraints to the start configuration and reset shells. - 'yes': Do not apply constraints to the start configuration and do not reset - shells, useful for exact coninuation and reruns. - """ shake_tol: float = 0.0001 """ Relative tolerance for SHAKE. Default 0.0001. """ lincs_order: int = 12 From 480ad7087930e12e58cca61ee008ec4adf2f242a Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Wed, 4 Sep 2024 16:48:30 +0200 Subject: [PATCH 22/23] Move settings that are irrelevant for EM to MDSimulationSettings --- .../protocols/gromacs_md/md_settings.py | 341 ++++++++++-------- 1 file changed, 187 insertions(+), 154 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 43da700..225a83d 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -43,14 +43,6 @@ class Config: set to that mMin. Default 1 (no mass scaling) """ - # # # Langevin dynamics # # # - ld_seed: int = -1 - """ - Integer used to initialize random generator for thermal noise for - stochastic and Brownian dynamics. When ld-seed is set to -1, - a pseudo random seed is used. Default -1. - """ - # # # Neighbor searching # # # cutoff_scheme: Literal["verlet"] = "verlet" """ @@ -154,115 +146,6 @@ class Config: Default 1e-5 """ - # # # Temperature coupling # # # - tcoupl: Literal[ - "no", "berendsen", "nose-hoover", "andersen", "andersen-massive", "v-rescale" - ] = "no" - """ - Temperature coupling options. Note that tcoupl will be ignored when the - 'sd' integrator is used. - Allowed values are: - 'no': No temperature coupling. - 'berendsen': Temperature coupling with a Berendsen thermostat to a bath - with temperature ref-t, with time constant tau-t. Several groups can be - coupled separately, these are specified in the tc-grps field separated - by spaces. This is a historical thermostat needed to be able to - reproduce previous simulations, but we strongly recommend not to use - it for new production runs. - 'nose-hoover': Temperature coupling using a Nose-Hoover extended ensemble. - The reference temperature and coupling groups are selected as above, - but in this case tau-t controls the period of the temperature - fluctuations at equilibrium, which is slightly different from a - relaxation time. For NVT simulations the conserved energy quantity is - written to the energy and log files. - 'andersen': Temperature coupling by randomizing a fraction of the particle - velocities at each timestep. Reference temperature and coupling groups - are selected as above. tau-t is the average time between randomization - of each molecule. Inhibits particle dynamics somewhat, but little or no - ergodicity issues. Currently only implemented with velocity Verlet, and - not implemented with constraints. - 'andersen-massive': Temperature coupling by randomizing velocities of all - particles at infrequent timesteps. Reference temperature and coupling - groups are selected as above. tau-t is the time between randomization - of all molecules. Inhibits particle dynamics somewhat, but little or - no ergodicity issues. Currently only implemented with velocity Verlet. - 'v-rescale': Temperature coupling using velocity rescaling with a - stochastic term (JCP 126, 014101). This thermostat is similar to - Berendsen coupling, with the same scaling using tau-t, but the - stochastic term ensures that a proper canonical ensemble is generated. - The random seed is set with ld-seed. This thermostat works correctly - even for tau-t =0. For NVT simulations the conserved energy quantity - is written to the energy and log file. - Default 'no' (for the default integrator, 'sd', this option is ignored). - """ - ref_t: FloatQuantity["kelvin"] = 298.15 * unit.kelvin - """ - Reference temperature for coupling (one for each group in tc-grps). - Default 298.15 * unit.kelvin - """ - - # # # Pressure coupling # # # - pcoupl: Literal["no", "berendsen", "C-rescale", "Parrinello-Rahman"] = "no" - """ - Options for pressure coupling (barostat). Allowed values are: - 'no': No pressure coupling. This means a fixed box size. - 'berendsen': Exponential relaxation pressure coupling with time constant - tau-p. The box is scaled every nstpcouple steps. This barostat does not - yield a correct thermodynamic ensemble; it is only included to be able - to reproduce previous runs, and we strongly recommend against using it - for new simulations. - 'C-rescale': Exponential relaxation pressure coupling with time constant - tau-p, including a stochastic term to enforce correct volume - fluctuations. The box is scaled every nstpcouple steps. It can be used - for both equilibration and production. - 'Parrinello-Rahman': Extended-ensemble pressure coupling where the box - vectors are subject to an equation of motion. The equation of motion - for the atoms is coupled to this. No instantaneous scaling takes place. - As for Nose-Hoover temperature coupling the time constant tau-p is the - period of pressure fluctuations at equilibrium. - Default 'no'. - """ - ref_p: FloatQuantity["bar"] = 1.01325 * unit.bar - """ - The reference pressure for coupling. The number of required values is - implied by pcoupltype. Default 1.01325 * unit.bar. - """ - refcoord_scaling: Literal["no", "all", "com"] = "no" - """ - Allowed values are: - 'no': The reference coordinates for position restraints are not modified. - 'all': The reference coordinates are scaled with the scaling matrix of the - pressure coupling. - 'com': Scale the center of mass of the reference coordinates with the - scaling matrix of the pressure coupling. The vectors of each reference - coordinate to the center of mass are not scaled. Only one COM is used, - even when there are multiple molecules with position restraints. - For calculating the COM of the reference coordinates in the starting - configuration, periodic boundary conditions are not taken into account. - Default 'no'. - """ - - # # # Velocity generation # # # - gen_vel: Literal["no", "yes"] = "yes" - """ - Velocity generation. Allowed values are: - 'no': Do not generate velocities. The velocities are set to zero when there - are no velocities in the input structure file. - 'yes': Generate velocities in gmx grompp according to a Maxwell - distribution at temperature gen-temp, with random seed gen-seed. - This is only meaningful with integrator=md. - Default 'yes'. - """ - gen_temp: FloatQuantity["kelvin"] = 298.15 * unit.kelvin - """ - Temperature for Maxwell distribution. Default 298.15 * unit.kelvin. - """ - gen_seed: int = -1 - """ - Used to initialize random generator for random velocities, when gen-seed is - set to -1, a pseudo random seed is used. Default -1. - """ - # # # Bonds # # # constraints: Literal[ "none", "h-bonds", "all-bonds", "h-angles", "all-angles" @@ -330,8 +213,6 @@ class Config: "rcoulomb", "rvdw", "ewald_rtol", - "ref_t", - "gen_temp", "shake_tol", "lincs_order", "lincs_iter", @@ -340,18 +221,17 @@ def must_be_positive_or_zero(cls, v): if v < 0: errmsg = ( "Settings nsteps, nstlist, rlist, rcoulomb, rvdw, ewald_rtol, " - "ref_t, gen_temp, shake_tol, lincs_order, and lincs_iter" + "shake_tol, lincs_order, and lincs_iter" f" must be zero or positive values, got {v}." ) raise ValueError(errmsg) return v - @validator("dt", "mass_repartition_factor") + @validator("mass_repartition_factor") def must_be_positive(cls, v): if v <= 0: errmsg = ( - "timestep dt, and mass_repartition_factor " - f"must be positive values, got {v}." + f"mass_repartition_factor must be positive values, got {v}." ) raise ValueError(errmsg) return v @@ -363,12 +243,6 @@ def must_be_between_3_12(cls, v): raise ValueError(errmsg) return v - @validator("dt") - def is_time(cls, v): - if not v.is_compatible_with(unit.picosecond): - raise ValueError("dt must be in time units " "(i.e. picoseconds)") - return v - @validator("rlist", "rcoulomb", "rvdw") def is_distance(cls, v): if not v.is_compatible_with(unit.nanometer): @@ -378,31 +252,6 @@ def is_distance(cls, v): ) return v - @validator("ref_t", "gen_temp") - def is_temperature(cls, v): - if not v.is_compatible_with(unit.kelvin): - raise ValueError( - "ref_t and gen_temp must be in temperature units" " (i.e. kelvin)" - ) - return v - - @validator("ref_p") - def is_pressure(cls, v): - if not v.is_compatible_with(unit.bar): - raise ValueError("ref_p must be in pressure units (i.e. bar)") - return v - - @validator("integrator") - def supported_integrator(cls, v): - supported = ["md", "sd", "steep"] - if v.lower() not in supported: - errmsg = ( - "Only the following sampler_method values are " - f"supported: {supported}, got {v}" - ) - raise ValueError(errmsg) - return v - class OutputSettings(SettingsBaseModel): """ " @@ -517,6 +366,17 @@ class EMSimulationSettings(SimulationSettings): 'steep': A steepest descent algorithm for energy minimization. """ + @validator("integrator") + def supported_integrator(cls, v): + supported = ["steep"] + if v.lower() not in supported: + errmsg = ( + "Only the following sampler_method values are " + f"supported: {supported}, got {v}" + ) + raise ValueError(errmsg) + return v + class EMOutputSettings(OutputSettings): """ @@ -544,6 +404,179 @@ class MDSimulationSettings(SimulationSettings): """ + # # # Langevin dynamics # # # + ld_seed: int = -1 + """ + Integer used to initialize random generator for thermal noise for + stochastic and Brownian dynamics. When ld-seed is set to -1, + a pseudo random seed is used. Default -1. + """ + + # # # Temperature coupling # # # + tcoupl: Literal[ + "no", "berendsen", "nose-hoover", "andersen", "andersen-massive", + "v-rescale" + ] = "no" + """ + Temperature coupling options. Note that tcoupl will be ignored when the + 'sd' integrator is used. + Allowed values are: + 'no': No temperature coupling. + 'berendsen': Temperature coupling with a Berendsen thermostat to a bath + with temperature ref-t, with time constant tau-t. Several groups can be + coupled separately, these are specified in the tc-grps field separated + by spaces. This is a historical thermostat needed to be able to + reproduce previous simulations, but we strongly recommend not to use + it for new production runs. + 'nose-hoover': Temperature coupling using a Nose-Hoover extended ensemble. + The reference temperature and coupling groups are selected as above, + but in this case tau-t controls the period of the temperature + fluctuations at equilibrium, which is slightly different from a + relaxation time. For NVT simulations the conserved energy quantity is + written to the energy and log files. + 'andersen': Temperature coupling by randomizing a fraction of the particle + velocities at each timestep. Reference temperature and coupling groups + are selected as above. tau-t is the average time between randomization + of each molecule. Inhibits particle dynamics somewhat, but little or no + ergodicity issues. Currently only implemented with velocity Verlet, and + not implemented with constraints. + 'andersen-massive': Temperature coupling by randomizing velocities of all + particles at infrequent timesteps. Reference temperature and coupling + groups are selected as above. tau-t is the time between randomization + of all molecules. Inhibits particle dynamics somewhat, but little or + no ergodicity issues. Currently only implemented with velocity Verlet. + 'v-rescale': Temperature coupling using velocity rescaling with a + stochastic term (JCP 126, 014101). This thermostat is similar to + Berendsen coupling, with the same scaling using tau-t, but the + stochastic term ensures that a proper canonical ensemble is generated. + The random seed is set with ld-seed. This thermostat works correctly + even for tau-t =0. For NVT simulations the conserved energy quantity + is written to the energy and log file. + Default 'no' (for the default integrator, 'sd', this option is ignored). + """ + ref_t: FloatQuantity["kelvin"] = 298.15 * unit.kelvin + """ + Reference temperature for coupling (one for each group in tc-grps). + Default 298.15 * unit.kelvin + """ + + # # # Pressure coupling # # # + pcoupl: Literal["no", "berendsen", "C-rescale", "Parrinello-Rahman"] = "no" + """ + Options for pressure coupling (barostat). Allowed values are: + 'no': No pressure coupling. This means a fixed box size. + 'berendsen': Exponential relaxation pressure coupling with time constant + tau-p. The box is scaled every nstpcouple steps. This barostat does not + yield a correct thermodynamic ensemble; it is only included to be able + to reproduce previous runs, and we strongly recommend against using it + for new simulations. + 'C-rescale': Exponential relaxation pressure coupling with time constant + tau-p, including a stochastic term to enforce correct volume + fluctuations. The box is scaled every nstpcouple steps. It can be used + for both equilibration and production. + 'Parrinello-Rahman': Extended-ensemble pressure coupling where the box + vectors are subject to an equation of motion. The equation of motion + for the atoms is coupled to this. No instantaneous scaling takes place. + As for Nose-Hoover temperature coupling the time constant tau-p is the + period of pressure fluctuations at equilibrium. + Default 'no'. + """ + ref_p: FloatQuantity["bar"] = 1.01325 * unit.bar + """ + The reference pressure for coupling. The number of required values is + implied by pcoupltype. Default 1.01325 * unit.bar. + """ + refcoord_scaling: Literal["no", "all", "com"] = "no" + """ + Allowed values are: + 'no': The reference coordinates for position restraints are not modified. + 'all': The reference coordinates are scaled with the scaling matrix of the + pressure coupling. + 'com': Scale the center of mass of the reference coordinates with the + scaling matrix of the pressure coupling. The vectors of each reference + coordinate to the center of mass are not scaled. Only one COM is used, + even when there are multiple molecules with position restraints. + For calculating the COM of the reference coordinates in the starting + configuration, periodic boundary conditions are not taken into account. + Default 'no'. + """ + + # # # Velocity generation # # # + gen_vel: Literal["no", "yes"] = "yes" + """ + Velocity generation. Allowed values are: + 'no': Do not generate velocities. The velocities are set to zero when there + are no velocities in the input structure file. + 'yes': Generate velocities in gmx grompp according to a Maxwell + distribution at temperature gen-temp, with random seed gen-seed. + This is only meaningful with integrator=md. + Default 'yes'. + """ + gen_temp: FloatQuantity["kelvin"] = 298.15 * unit.kelvin + """ + Temperature for Maxwell distribution. Default 298.15 * unit.kelvin. + """ + gen_seed: int = -1 + """ + Used to initialize random generator for random velocities, when gen-seed is + set to -1, a pseudo random seed is used. Default -1. + """ + + @validator( + "ref_t", + "gen_temp", + ) + def must_be_positive_or_zero(cls, v): + if v < 0: + errmsg = ( + "Settings ref_t, and gen_temp must be zero or positive values," + f" got {v}." + ) + raise ValueError(errmsg) + return v + + @validator("dt") + def must_be_positive(cls, v): + if v <= 0: + errmsg = ( + f"timestep dt must be positive values, got {v}." + ) + raise ValueError(errmsg) + return v + + @validator("dt") + def is_time(cls, v): + if not v.is_compatible_with(unit.picosecond): + raise ValueError("dt must be in time units " "(i.e. picoseconds)") + return v + + @validator("ref_t", "gen_temp") + def is_temperature(cls, v): + if not v.is_compatible_with(unit.kelvin): + raise ValueError( + "ref_t and gen_temp must be in temperature units" " (i.e. " + "kelvin)" + ) + return v + + @validator("ref_p") + def is_pressure(cls, v): + if not v.is_compatible_with(unit.bar): + raise ValueError("ref_p must be in pressure units (i.e. bar)") + return v + + @validator("integrator") + def supported_integrator(cls, v): + supported = ["md", "sd"] + if v.lower() not in supported: + errmsg = ( + "Only the following sampler_method values are " + f"supported: {supported}, got {v}" + ) + raise ValueError(errmsg) + return v + + class NVTSimulationSettings(MDSimulationSettings): """ Settings for MD simulation in the NVT ensemble. From e01fdd87d13f9a5cc7afd973919e455f50490927 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 4 Sep 2024 14:51:14 +0000 Subject: [PATCH 23/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../protocols/gromacs_md/md_settings.py | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/openfe_gromacs/protocols/gromacs_md/md_settings.py b/openfe_gromacs/protocols/gromacs_md/md_settings.py index 225a83d..ce04d1b 100644 --- a/openfe_gromacs/protocols/gromacs_md/md_settings.py +++ b/openfe_gromacs/protocols/gromacs_md/md_settings.py @@ -230,9 +230,7 @@ def must_be_positive_or_zero(cls, v): @validator("mass_repartition_factor") def must_be_positive(cls, v): if v <= 0: - errmsg = ( - f"mass_repartition_factor must be positive values, got {v}." - ) + errmsg = f"mass_repartition_factor must be positive values, got {v}." raise ValueError(errmsg) return v @@ -403,7 +401,6 @@ class MDSimulationSettings(SimulationSettings): Default 0.002 * unit.picosecond """ - # # # Langevin dynamics # # # ld_seed: int = -1 """ @@ -414,8 +411,7 @@ class MDSimulationSettings(SimulationSettings): # # # Temperature coupling # # # tcoupl: Literal[ - "no", "berendsen", "nose-hoover", "andersen", "andersen-massive", - "v-rescale" + "no", "berendsen", "nose-hoover", "andersen", "andersen-massive", "v-rescale" ] = "no" """ Temperature coupling options. Note that tcoupl will be ignored when the @@ -538,9 +534,7 @@ def must_be_positive_or_zero(cls, v): @validator("dt") def must_be_positive(cls, v): if v <= 0: - errmsg = ( - f"timestep dt must be positive values, got {v}." - ) + errmsg = f"timestep dt must be positive values, got {v}." raise ValueError(errmsg) return v @@ -554,8 +548,7 @@ def is_time(cls, v): def is_temperature(cls, v): if not v.is_compatible_with(unit.kelvin): raise ValueError( - "ref_t and gen_temp must be in temperature units" " (i.e. " - "kelvin)" + "ref_t and gen_temp must be in temperature units" " (i.e. " "kelvin)" ) return v