|
apoCHARMM 1.0.0
High-performance molecular dynamics simulations on GPUs
|
|
|
apoCHARMM 1.0.0
High-performance molecular dynamics simulations on GPUs
|
|
The CUDA-integrator subsystem advances a CharmmContext through molecular-dynamics steps on a CUDA device. CudaIntegrator provides the common context, time-step, restart, subscriber, step-counter, neighbor-list, center-of-mass, and CUDA-stream lifecycle. It is a polymorphic base that delegates initialization and each dynamics step to a concrete integrator.
The currently documented concrete choices are:
The native base class can be constructed directly, but its default initialization, restart, descriptor, and single-step hooks report ApoCharmmErrorCode::NotImplemented. Normal user code constructs one of the concrete classes.
The following Python example propagates a small system with the stochastic Langevin thermostat:
The time step is 0.001 ps, the target temperature is 300 K, and the thermostat friction is 1 ps^-1. After successful propagation, the final line prints the finite running-average temperature maintained by the integrator.
The corresponding native workflow is:
Repository-relative test/data paths are suitable for repository tests and examples. Applications should supply their own topology, parameter, and coordinate paths.
Public time-step arguments use picoseconds and must be finite and greater than zero. The base implementation stores the value internally in AKMA time by dividing by 0.0488882129; getTimeStep() converts it back to picoseconds.
Configure the time step and concrete-integrator parameters before attaching the context. setCharmmContext() may be called only once. It retains the context, resizes the base per-atom work arrays, constructs any configured holonomic-constraint helper, and immediately invokes the concrete initialization hook. A failure clears the retained context and constraint helper, but previously resized work storage can remain observable.
A context used for propagation must have a non-null, initialized ForceManager. The context must also contain the coordinates, velocity/inverse-mass state, atom count, degree-of-freedom count, box state, and other force-manager inputs required by the selected force configuration. propagate() rejects a non-positive number of steps and checks context and force-manager readiness before changing the trajectory.
For CudaLangevinPistonIntegrator, call setCrystalType() before attaching the context. The supported layouts are:
CRYSTAL::CUBIC: one isotropic box degree of freedom.CRYSTAL::TETRAGONAL: one shared X/Y degree of freedom and one Z degree of freedom.CRYSTAL::ORTHORHOMBIC: independent X, Y, and Z degrees of freedom.Setting the crystal type allocates and resets all crystal-dependent piston arrays. A custom Langevin-piston mass must contain exactly one, two, or three values according to that selected layout.
Attaching a context initializes the concrete integrator automatically. initialize() repeats concrete initialization and can overwrite derived displacement, random-number, kinetic, or extended-variable state. Use initializeFromRestartFile() after attaching the context when restoring a saved trajectory.
setTimeStep() changes only the base time-step scalar. It does not rebuild an existing holonomic-constraint helper or recalculate every derived coefficient and displacement history. Treat the time step as immutable after context attachment.
A native integrator retains its CharmmContext with a std::shared_ptr. It also retains every subscribed Subscriber with a std::shared_ptr. Subscription installs context and integrator backlinks in the subscriber, so the native graph can contain a strong-reference cycle. Unsubscribe objects explicitly before releasing the last external integrator and subscriber owners.
Each native integrator owns two CUDA stream handles through shared pointers. The concrete stochastic integrators additionally own raw CUDA allocations for Philox random-number states. Their destructors are non-throwing and suppress CUDA cleanup errors.
The C ABI concrete handle owns a shared native concrete integrator and embeds an apo_cuda_integrator base view. Conversion functions ending in _as_cuda_integrator return a borrowed pointer to that embedded view. The borrowed view has no destroy function and becomes invalid when its concrete handle is destroyed.
apo_cuda_integrator_set_charmm_context() borrows the C context handle for one call but retains the shared native context on success. The source context handle may then be destroyed without invalidating the integrator. Subscription similarly retains the native subscriber, but the concrete C subscriber handle must remain alive when its embedded base view will later be passed to unsubscription.
Each Python concrete wrapper owns its concrete C handle and stores a borrowed base-integrator view. The base wrapper retains the supplied Python context and successfully subscribed Python subscriber wrappers. close() destroys the concrete C handle, invalidates both handle properties, and clears those Python references. It does not implicitly unsubscribe native subscribers; unsubscribe before closure.
Native, C ABI, and Python objects provide no locking. Do not overlap mutation, initialization, propagation, restart loading, subscription, unsubscription, or destruction involving the same integrator or context from multiple host threads.
The public quantities whose units are established by the current implementation are:
| Quantity | Public representation | Unit |
|---|---|---|
| Time step | Scalar | picoseconds |
| Internal time step | Native protected scalar | AKMA time |
| Reference and measured temperature | Scalar | kelvin |
| Thermostat and piston friction | Scalar | inverse picoseconds |
| Kinetic energy | Two-element native container | kilocalories per mole |
| Pressure tensor | Nine row-major values | atmospheres |
| Pressure scalar | One value | atmospheres |
| Box and piston position state | Scalar arrays | angstroms |
| Piston on-step velocity | Scalar array | angstroms per AKMA time |
| Piston delta pressure | Scalar array | atmospheres per angstrom |
| Instantaneous surface tension | Scalar | dyne per centimeter |
| Seeds, sequence positions, frequencies, and iteration counts | Integer | dimensionless |
The exact dimensional convention for the public Nose-Hoover coupling mass, Nose-Hoover extended-variable velocity/increment values, and Langevin-piston mass is not stated by the current repository. Their setters and getters therefore document only validation, shape, copying, and the native numeric convention.
The base displacement arrays are CudaContainer objects with one double4 record per atom. The x, y, and z components are angstrom-valued displacement or reference-coordinate components. The integrators do not assign public meaning to the w component. Native accessors return borrowed container aliases and perform no implicit host/device transfer.
Nose-Hoover and Langevin-piston kinetic-energy and running-temperature arrays have two elements:
The native instantaneous-temperature methods select element 1 when useOldTemperature(true) is active and element 0 otherwise. The current Nose-Hoover and Langevin-piston C ABI running-average getters select the opposite indices: they return element 0 when the old-temperature flag is true and element 1 when it is false. Their Python wrappers inherit that current C ABI behavior.
The Langevin thermostat always reports element 0 through its C ABI and Python running-average and instantaneous-temperature methods.
Piston pressure tensors are flat, row-major 3-by-3 arrays in this order:
The Python piston wrapper accepts an outer iterable of inner iterables and flattens all values in iteration order. Native validation requires exactly nine total values but the current wrapper does not independently require three rows of three values. Tensor getters return newly allocated three-by-three Python lists.
The Langevin thermostat maintains one Philox state per atom. The Langevin piston maintains one state per active crystal degree of freedom. Native serialized RNG-state strings are opaque implementation representations; callers may save and restore them through the native API but should not parse or modify their contents.
Verified validation and lifecycle failures use ApoCharmmError:
ApoCharmmErrorCode::InvalidArgument reports invalid scalar ranges, null or duplicate attachments, invalid array lengths, unsupported crystal values, incompatible restart state, and invalid propagation counts.ApoCharmmErrorCode::NotInitialized reports missing context, force-manager, crystal, RNG, or other required state.ApoCharmmErrorCode::NotImplemented reports direct use of unimplemented base hooks and unsupported subscriber operations.ApoCharmmErrorCode::Runtime reports restart-file open, format, parse, and state-shape failures, NaN detection, and other explicitly categorized runtime failures.ApoCharmmErrorCode::Cuda reports checked CUDA allocation, copy, stream, synchronization, and kernel-launch failures.Host allocation and string construction can additionally propagate standard exceptions such as std::bad_alloc, std::length_error, or std::system_error. Subscriber implementations retain a small number of legacy std::invalid_argument and std::runtime_error paths.
Initialization, restart loading, propagation, and multi-subscriber operations are not transactional. State changes, completed steps, updated counters, and written output from before a failure remain observable.
Every status-returning integrator function clears the calling thread's previous diagnostic at entry. Native ApoCharmmErrorCode values map to their corresponding apo_status; std::invalid_argument maps to APO_STATUS_INVALID_ARGUMENT; other C++ exceptions map to APO_STATUS_RUNTIME_ERROR.
After a nonzero status, call apo_last_error() immediately on the same thread. The returned diagnostic pointer is borrowed and remains valid only until the next diagnostic-changing C ABI call on that thread or thread exit.
Owned-handle output slots are set to NULL before later construction or validation. Scalar and fixed-buffer getters write their output only after handle, pointer, shape, and transfer operations succeed unless the individual function documents a different count-output rule. Destroy functions accept NULL, do not throw, and preserve an existing diagnostic when destruction succeeds.
Python wrapper checks raise TypeError for rejected object types or native scalar conversions and ValueError for explicit integer-range checks. Access through a closed wrapper raises RuntimeError. File-path encoding can raise TypeError, ValueError, or OSError according to os.fsencode().
Every nonzero C ABI status is raised as apocharmm.ApoCharmmError with the native status and diagnostic. Python array inputs are copied into temporary C buffers. Tensor, mass, and scalar outputs are returned as new Python values and do not alias native storage.
propagate() resets the context neighbor list before the loop, sets the per-call neighbor-list counter to zero, stores the requested step count, and adds the complete request to the total-step counter before executing the first step. A later failure therefore leaves the total requested-step count greater than the number of completed steps.
The current propagated-step counter is incremented after each successful concrete step and wraps from INT_MAX to zero. Calling propagateOneStep() directly does not increment these counters, run scheduled subscribers, reset the neighbor list, or perform the periodic NaN check.
Subscriber reporting intervals are copied into a parallel scheduling array at subscription time. Changing a subscriber's own reporting frequency afterward does not change the cached schedule. Each call to propagate() numbers its local subscriber schedule from one.
The default nonbonded-list update interval is 20 steps and the default center-of-mass-removal interval is 1000 steps. Because the current-step counter starts at zero, each concrete integrator performs its center-of-mass operation during the first propagated step. These operations can synchronize the integrator stream and transfer velocity or displacement storage between device and host.
Concrete initialization and each complete concrete step synchronize the private integrator stream at their documented completion boundary. Intermediate work is issued asynchronously on that stream. Context force calculations and explicit CudaContainer transfers can introduce additional stream or device-wide synchronization. The separately allocated integrator memcpy stream is not used by the three inspected concrete implementations.
The Langevin thermostat advances its recorded RNG sequence position by four per propagated step because each per-atom Philox call consumes a curand_normal4 group. The Langevin piston advances its position once per predictor-corrector iteration, so one propagated step advances it by the configured maximum predictor-corrector iteration count.
Changing a Langevin thermostat seed or sequence position reinitializes its per-atom RNG states when a context is attached. Changing a piston seed or sequence position reinitializes its RNG states only after a crystal type has established the piston degree-of-freedom count.
Mutable native container accessors bypass shape, coherence, and initialization invariants. Resizing or replacing kinetic-energy, temperature, pressure, extended-variable, or displacement arrays can cause later kernels or C ABI getters to read invalid layouts. They are contributor-level escape hatches, not ordinary configuration APIs.
The current implementation has several material limitations:
getAverageOldTemperature() piston declarations have no definitions.The public native declarations are in:
Their implementations and CUDA kernels are in the correspondingly named files under src/. The public C ABI is declared under capi/include/apocharmm_c/, implemented under capi/src/apocharmm_c/, and backed by private handle structs under capi/private_include/apocharmm_c/detail/. Python wrappers under python/apocharmm/ use ctypes and those C ABI entry points.
The native lifecycle is:
A concrete extension derives from CudaIntegrator and implements initializeImpl(), initializeFromRestartFileImpl(), and propagateOneStepImpl(). The base final methods enforce the common lifecycle and readiness boundary. An extension must preserve the invariant that m_CoordsRef, m_CoordsDelta, and m_CoordsDeltaPrevious contain one double4 record per context atom after successful attachment.
The base owns the scheduling arrays. Entry i of m_Subscribers must always correspond to entry i of m_ReportFreqList. The mutable native accessors to those arrays are retained for compatibility but can violate this invariant.
The deterministic Nose-Hoover and piston paths derive two kinetic estimators from previous-half-step displacement, on-step velocity, and next-half-step displacement. A predictor-corrector loop updates the Nose-Hoover extended variable. The piston path additionally computes kinetic and virial pressure, updates crystal degrees of freedom, scales coordinates and box dimensions, and repeats the piston/pressure correction for the configured iteration count.
The stochastic thermostat updates positions with Philox normal variates, applies holonomic constraints when enabled, derives velocities from corrected displacements, and computes both kinetic estimators. The stochastic piston uses one Philox state per crystal degree of freedom and applies random forcing inside the predictor-corrector loop.
Host/device flow is explicit. Kernels consume context and integrator device arrays. C ABI scalar and buffer getters generally call transferToHost() before reading a native CudaContainer. Restart readers parse host strings and arrays, then transfer restored state to the device. No automatic dirty-state tracking exists between mutable native aliases and later consumers.
The native ApoCharmmError boundary is inside the C++ implementation. The C ABI guard translates all escaping C++ failures into apo_status and thread-local diagnostics. The Python shared-library wrapper installs a status error checker, so nonzero statuses become ApoCharmmError before wrapper methods continue to their post-call ownership updates.
Performance-sensitive operations include context attachment, holonomic setup, RNG allocation and initialization, force calculation, neighbor-list rebuilding, full per-atom device copies, pressure reductions, predictor-corrector loops, subscriber transfers, and stream synchronization. Repeated parameter changes that trigger reinitialization should not be placed in a propagation loop.
The focused native regression files are:
Direct C ABI coverage exists in:
There is currently no direct C ABI Nose-Hoover integrator test. Python API coverage is in the four python_api_cuda_*integrator.py files under test/pytest/. End-to-end examples are:
1.9.8